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ABSTRACT 

More  accurate,  more  reliable  and  more  efficient  methods  to  calculate  the  bistatic  scattering  of 
electromagnetic  fields  from  absorbing  dielectric  objects  have  been  developed.  While  the  solution  of  three 
dimensional  scattering  problems  has  not  been  achieved,  the  methods  developed  here  provide  for  the  first  time 
finite  element  solutions  of  closed  three-dimensional  electromagnetic  fields  free  of  the  spurious  solutions  that 
have  plagued  previous  procedures.  These  stable  finite  elements  are  defined  in  three  varieties:  (1)  mixed-order 
rectangular  parallelepiped  elements,  (2)  edge-based  tangential  vector  finite  elements,  and  (3)  derivative 
continuous  C1  elements.  Each  are  shown  to  give  only  physically  correct  solutions.  In  addition,  a  new 
procedure  called  the  transfinite  element  method  is  presented  for  the  solution  of  open  or  unbounded 
electromagnetic  scattering  problems.  This  method  provides  a  one  step  procedure  to  compute  scattered  fields  in 
electromagnetics.  Prior  to  the  transfinite  element  method,  scattering  problems  required  that  many  time  steps 
or  many  modes  be  used  to  produce  each  field  solution.  The  transfinite  element  method  is  applied  here  only  in 
two  dimensions.  Finally,  a  new  high  efficiency  algorithm  is  proposed  to  compute  electromagnetic  scattering 
over  a  specified  frequency  range.  In  this  algorithm,  frequencies  are  selected  adaptively  to  ensure  a  given 
accuracy  in  minimum  computation  time  for  the  specified  frequency  range.  The  individual  procedures  are 
illustrated  with  a  variety  of  examples.  Thus,  a  great  deal  of  progress  in  the  solution  of  electromagnetic 
scattering  problems  has  been  achieved,  although  the  major  program  goal  of  solving  for  the  electromagnetic 
field':  scattered  by  large  three  dimensional  objects  has  not  been  accomplished. 
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SECTION  1  INTRODUCTION 

A  program  of  work  has  been  completed  to  develop  more  accurate,  more  reliable  and  more  efficient 
numerical  procedures  for  the  computation  of  electromagnetic  fields  scattered  by  absorbing  dielectric  objects. 
This  work  has  focused  on  the  application  of  the  finite  element  method  to  this  problem  since  the  finite  element 
method  is  well  suited  to  modeling  problems  involving  complex  shapes  and  material  properties.  It  involved  the 
solution  of  two  separate  fundamental  problems  in  electromagnetics:  (1)  How  does  one  approximate 
electromagnetic  fields  over  finite  element  domains  such  that  stable  solutions  are  produced?  and  (2)  How  does 
one  interface  the  open  boundary  condition  encountered  in  scattering  problems  to  the  finite  element  solution  in 
an  accurate  and  efficient  manner? 

The  first  application  of  the  finite  element  method  in  electromagnetics  was  published  by  Silvester  in 
1969[1].  In  this  application,  the  finite  element  method  was  used  to  solve  for  the  resonant  electromagnetic 
fields  in  homogeneous  waveguides.  Subsequently,  Silvester  and  his  coworkers  published  a  series  of  papers 
[2-6]  extending  the  realm  of  application  of  the  finite  element  method  in  high  frequency  electromagnetics. 
However,  it  soon  became  apparent  that  the  category  of  problem  with  which  the  finite  element  method  was 
stable  was  quite  limited:  the  finite  element  method  only  gave  good  answers  to  problems  in  which  a  single 
scalar  variable  was  the  unknown.  Three-dimensional  problems  or  two-dimensional  problems  that  required  the 
solution  of  vector  field  quantities  gave  spurious  solutions[3-6].  In  fact,  the  subject  of  spurious  modes  in 
finite  element  solutions  of  electromagnetics  problems  has  generated  considerable  literature  in  recent 
years[7-ll]. 

Parallel  to  the  work  of  Silvester  and  others  on  the  application  of  the  finite  element  method  to  closed 
high-frequency  electromagnetics  problems,  Mei  and  his  coworkers  developed  a  procedure  called  the 
unimoment  method  for  solving  scalar  variable  unbounded  electromagnetic  field  problems  using  the  finite 
element  method[12-13].  This  procedure  makes  it  possible  to  solve  complicated  scattering  problems  by  placing 
the  scatterer  inside  a  circular  boundary,  discretizing  the  region  inside  the  circle  with  finite  elements,  and 
coupling  the  finite  element  solution  to  the  analytical  solution  outside  the  circle  through  interface  conditions. 
Unfortunately,  however,  this  coupling  requires  that  a  full  set  of  modal  finite  element  solutions  be  evaluated 
such  that  each  finite  element  solution  and  its  derivative  match  the  corresponding  free  space  solution  on  the 
circular  interface.  This  process  is  relatively  time  consuming  because  many  finite  element  problems  need  to  be 
solved  to  obtain  a  single  scattered  field  solution  and  it  is  relatively  inaccurate  because  numerical  derivatives  are 
employed  in  the  formulation. 

The  goal  in  this  work  is  to  develop  numerical  procedures  by  which  relatively  large  three 
dimensional  dielectric  scatterers  can  be  modeled.  To  achieve  this  goal,  two  new  developments  were  required: 
First,  it  was  essential  to  develop  numerical  procedures  for  modeling  three  dimensional  electromagnetic  field 
problems  in  a  stable  way.  Second,  more  efficient  numerical  procedures  fen’  setting  up  and  solving  open 
electromagnetic  field  problems  were  needed 
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This  report  documents  progress  we  have  made  towards  achieving  these  goals.  With  respect  to  the 
problem  of  spurious  solutions  in  three  dimensional  electromagnetic  field  calculations,  we  have  developed  three 
new  finite  element  procedures  by  which  the  problem  of  spurious  solutions  in  electromagnetics  is  eliminated. 
These  procedures  are:  (1)  mixed-order  finite  elements  over  rectangular  parallelepipeds,  (2)  tangentially 
continuous  vector  finite  elements  over  tetrahedrons,  and  (3)  C*  derivative  continuous  finite  elements  over 
tetrahedrons.  With  respect  to  the  problem  of  more  efficient  modeling  of  open  electromagnetics  problems,  we 
have  developed  a  new  procedure  called  the  transfmite  element  method  that  is  orders  of  magnitude  more  efficient 
and  more  accurate  than  the  unimoment  method.  We  have  also  developed  an  adaptive  spectral  modeling 
procedure  that  speeds  up  the  solution  of  scattering  problems  over  a  range  of  frequencies  by  at  least  an  order  of 
magnitude.  Although  time  did  not  permit  us  to  put  everything  together  to  solve  three  dimensional  scattering 
problems,  we  have  developed  the  major  theoretical  components  required  for  this  purpose. 


1.1  MOTIVATION 


1.1.1  INTEGRAL  VERSUS  DIFFERENTIAL  METHODS 

Most  common  methods  for  solving  electromagnetic  scattering  problems  are  based  on  integral 
equation  formulations.  In  general,  these  equations  are  discretized  by  using  pulse  or  linear  piecewise 
polynomial  basis  functions  and  are  solved  by  using  Galerkin's  method.  Originated  by  Roger  F.  Harrington, 
and  dubbed  the  method  of  moments,  such  procedures  for  solving  the  electromagnetic  field  equations  have  been 
applied  in  both  two  and  three  dimensional  scattering.  Although  these  procedures  have  been  successful  in 
solving  two  dimensional  scattering  problems,  their  success  in  modeling  three  dimensional  scattering  has  been 
more  limited. 

An  alternative  to  solving  electromagnetic  field  problems  using  integral  equation  methods  is  to  solve 
these  problems  by  using  methods  based  on  differential  formulations  of  the  electromagnetic  field.  The  two  best 
known  methods  to  solve  differential  equations  are  the  finite  difference  and  the  finite  element  method.  Although 
there  are  some  advantages  to  differential  methods  •  such  as  much  simpler  expressions  of  the  field  equations  - 
these  methods  have  been  seldom  used  to  solve  electromagnetic  scattering  problems.  Indeed,  until  recently,  in 
the  United  States  there  were  only  two  research  groups  working  on  modeling  electromagnetic  scattering  using 
differential  methods:  (K.  K.  Mei  at  the  University  of  California  at  Berkeley  and  Alan  Taflove  at  Northwestern 
University).  There  are  two  obvious  reasons  for  this:  (1)  the  number  of  data  points  required  by  differential 
methods  appears  at  first  glance  to  be  daunting,  and  (2)  representing  open,  unbounded  field  problems  is  difficult 
with  differential  methods.  With  integral  methods  these  issues  are  less  important:  only  the  surfaces  of  objects 
need  to  be  discretized  with  integral  methods  so  that  relatively  few  data  points  are  required.  Further,  the  far  field 
boundary  conditions  are  embedded  in  the  Green's  functions  used  in  integral  formulations  so  that  open 
problems  are  easily  addressed. 
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While  it  thus  appears  that  integral  methods  have  the  advantage  in  solving  electromagnetic  scattering 
problems,  this  is  in  fact  not  the  case.  As  described  in  more  detail  below,  a  large  number  of  methods  to  solve 
open  field  problems  by  truncating  the  finite  difference  or  finite  element  mesh  have  been  developed.  And, 
surprisingly,  the  advantage  of  less  computer  storage  requirements  and  of  less  computer  time  requirements  is 
actually  on  the  side  of  differential  methods.  This  is  surprising  because  differential  methods  require  many  more 
node  points  dm  integral  methods.  However,  all  of  the  interactions  in  a  differential  method  are  local  in  the 
sense  that  all  information  in  a  differential  model  is  passed  directly  from  one  element  to  the  next.  The 
differential  system  is  tied  together  by  these  local  interactions,  and  the  response  of  individual  elements  is 
obtained  by  solving  a  matrix  equation.  With  integral  methods,  information  may  pass  directly  from  one  side  of 
the  structure  to  the  other.  The  entire  integral  structure  is  tied  together  analytically  from  the  very  beginning  by 
evaluating  the  interaction  that  occurs  directly  between  any  two  points.  This  process  is  not  only  analytically 
challenging,  it  is  also  computationally  intensive. 

Thus,  although  differential  methods  are  at  a  disadvantage  because  the  full  space  must  be  modeled  and 
not  just  the  object  surfaces,  they  do  have  the  advantage  that  global  interactions  are  obtained  numerically  by 
means  of  the  matrix  solution  process  rather  that  by  using  analytical  expressions.  Since  computing  the 
interactions  between  different  surface  elements  with  integral  methods  is  expensive,  integral  methods  are  at  a 
disadvantage  in  this  respect  compared  to  differential  methods.  Further,  each  element  in  an  integral  method 
interacts  with  many  other  elements  so  that  the  matrix  equation  to  be  solved  to  compute  scattering  is  full.  In 
contrast,  the  matrices  arising  in  differential  methods  are  sparse  since  each  element  interacts  only  with  its 
immediate  neighbors. 

The  efficiency  of  numerical  methods  in  solving  electromagnetic  scattering  problems  is  determined 
by  several  factors:  (1)  the  efficiency  of  the  mesh  used  to  represent  the  objects  and  the  solution  variables,  (2) 
the  difficulty  and  cost  of  mesh  generation,  (3)  the  cost  of  generating  the  solution  matrix,  (4)  the  matrix 
storage  and  solution  time  requirements  of  the  algorithm,  and  (5)  the  robustness  of  computer  implementations 
with  respect  to  geometric  and  material  complexity.  These  factors  are  examined  next: 

1 .  Mesh  efficiency.  An  important  consideration  in  evaluating  the  efficiency  of  both  differential  and 
integral  solution  methods  is  the  type  of  mesh  required  for  its  operation.  The  application  of  some 
numerical  methods  is  restricted  to  regular  or  rectangular  meshes  while  other  methods  allow  irregular 
meshes  with  arbitrary  orientations  and  sizes.  Except  in  special  cases,  a  mesh  that  conforms  to  the 
boundaries  of  the  objects  that  are  modeled  and  that  is  refined  locally  in  regions  of  rapid  field 
variation  is  much  more  efficient  in  modeling  electromagnetic  field  problems  than  is  a  regular, 
rectangular  mesh.  The  gain  in  efficiency  depends  on  a  number  of  factors  such  as  the  problem 
geometry  and  the  method  used  to  construct  the  mesh  and  is  therefore  difficult  to  quantify. 
Fortunately,  however,  methods  for  constructing  optimal  triangular  meshes  in  two  dimensions  and 
optimal  tetrahedral  meshes  in  three  dimensions  have  been  developed  [29-33].  These  methods  are 
based  on  the  concept  of  Delaunay  tessellation  that  guarantees  that  the  mesh  employed  is  optimal 
with  respect  to  the  element  angles  regardless  of  the  point  spacing.  Furthermore,  with  irregular 
meshes,  it  is  possible  to  construct  the  mesh  adaptively  [29-32].  In  adaptive  mesh  generation, 
elements  are  refined  iteratively  in  the  areas  having  the  largest  error  until  all  elements  in  the  mesh 
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have  approximately  the  same  error.  In  terms  of  accuracy  for  a  given  mesh  size,  adaptive  Delaunay 
mesh  generation  is  the  most  efficient  procedure,  followed  by  non-adaptive  Delaunay  mesh 
generation,  other  irregular  mesh  generation  algorithms,  and,  last,  regular  mesh  generation.  The 
advantage  of  regular  mesh  generation  is,  of  course,  that  it  is  very  easy  and  inexpensive  to  compute. 
We  believe,  however,  that  the  modeling  efficiency  obtained  by  using  irregular  meshes  outweighs 
the  increased  algorithmic  complexity  of  these  procedures  in  most  cases. 

2.  Difficulty  and  cost  of  mesh  generation.  Since  only  the  surfaces  of  objects  need  to  be 
discretized  with  integral  methods  while  the  entire  volume  must  be  discretized  with  differential 
methods,  it  appears  at  first  that  it  is  easier  to  create  a  mesh  with  integral  methods.  However,  it 
turns  out  that  this  is  not  true  if  the  goal  is  to  obtain  irregular  meshes  for  arbitrary  objects  having 
complicated  shapes.  (Creating  regular  meshes  corresponding  to  simple  shapes  is  easy  with  both 
differential  and  integral  methods).  The  reason  for  this  is  that  the  Delaunay  approach  to  mesh 
generation  is  well  known  and  can  be  automated  in  both  two  and  three  dimensions.  No  direct, 
automatic,  optimal  procedure  to  generate  the  surface  mesh  of  arbitrary  objects  is  known.  At  this 
time,  the  only  known  way  to  create  an  optimal  surface  mesh  for  an  arbitrary  object  is  to  first  create 
the  three  dimensional  volume  Delaunay  mesh  for  the  object  and  then  to  extract  the  corresponding 
surface  mesh.  Thus,  the  difficulty  and  cost  of  mesh  generation  is  roughly  the  same  for  both 
differential  and  integral  methods. 

3 .  Matrix  set  up  cost.  Differential  methods  have  a  very  great  advantage  over  integral  methods 
with  respect  to  matrix  evaluation  cost.  When  they  are  available,  the  formulas  to  compute  matrix 
entries  with  differential  methods  are  simple  and  just  involve  arithmetic  expressions.  In  comparison, 
relatively  complicated  analytic  expressions  are  required  to  set  up  matrices  with  integral  methods. 
The  net  result  is  that  matrix  evaluation  cost  is  negligible  with  differential  methods,  but  is  often 
substantial  with  integral  methods. 

4.  Matrix  storage  and  solution  time  requirements.  The  analysis  of  matrix  storage  and 
solution  time  requirements  is  complicated  because  several  different  solution  procedures  exist  In  the 
past  decade,  the  use  of  the  preconditioned  conjugate  gradient  algorithm  has  had  a  profound  effect  on 
both  differential  and  integral  methods.  The  advantage  of  the  preconditioned  conjugate  gradient 
algorithm  is  that  with  differential  methods  matrix  solution  time  increases  linearly  with  matrix  size, 
compared  to  order  solution  times  for  solving  full  matrix  equations.  In  practice  with  non-adaptive 
differential  meshing  procedures,  solution  time  grows  as  with  both  two  and  three  dimensional 
geometries  [28].  However,  if  adaptive  finite  element  mesh  generation  procedures  are  employed,  the 
theoretical  ideal  limit  is  achieved  [29].  Thus,  in  the  finite  element  discretization  of  a  three 
dimensional  problem  with  an  n  X  n  X  n  mesh,  the  matrix  size  N  grows  as  N  *  n^  and  storage  and 
solution  time  for  these  sparse  matrices  increases  linearly.  The  net  result  is  that  computer  storage 
requirements  and  solution  times  grow  as  N  =  with  differential  methods.  This  value  can  be 
achieved  in  practice  by  using  adaptive  or  multigrid  methods.  (The  primary  difference  between 
multigrid  and  adaptive  methods  are  that  multigrid  methods  employ  rectangular  meshes  while 
adaptive  methods  allow  irregular  meshes;  they  are  similar  in  that  several  levels  of  meshes  are 


required).  If  non-adaptive  mesh  generation  procedures  are  employed,  then  the  performance  of  the 
algorithm  with  respect  to  solution  time  deteriorates  to  approximately  A 

With  boundary  element  methods,  only  the  elements  on  the  surface  of  three  dimensional  objects  need 
to  be  modeled  and  the  corresponding  mesh  grows  as  N  =  n^.  The  straightforward  application  of 
direct  solvers  to  the  corresponding  full  matrices  results  in  storage  and  solution  times  growing  as 
N^.  Thus,  with  integral  methods  and  full  matrix  solution  algorithms,  computer  requirements  grow 
as  n^.  The  application  of  the  conjugate  gradient  algorithm  to  these  matrices  lowers  the  solution 
time  requirement  to  n^  and  recent  work  by  Catedra,  Gago  and  Nuno  [34]  shows  that  combining  the 
conjugate  gradient  method  and  the  fast  Fourier  transform  reduces  this  further  to  n 2-667  i0g  n  xhis 
last  result  however  only  applies  to  problems  having  a  uniform  mesh  and  hence  applies  only  to 
problems  having  simple  shapes.  As  stated  earlier,  with  complicated  geometries  it  is  very  inefficient 
to  use  a  uniform  mesh  and  hence  the  fast  Fourier  transform  conjugate  gradient  method  is  not 
competitive  in  these  cases. 

5 .  Robustness.  An  important  consideration  in  evaluating  numerical  algorithms  for  modeling 
electromagnetic  scattering  is  robustness  with  respect  to  the  complexity  of  the  geometry  and  of  the 
material  properties  of  the  scatterers.  Here  differential  methods  have  a  clear  advantage:  provided  that 
the  same  mesh  is  used  in  both  cases,  with  differential  methods  it  makes  relatively  little  difference 
whether  the  material  characteristic  of  the  entire  problem  domain  is  homogeneous  or  whether  the 
material  characteristic  of  every  single  element  in  the  mesh  is  different  Thus,  although  differential 
methods  solution  requirements  grow  as  the  size  of  an  object  grows,  they  are  relatively  insensitive  to 
the  complexity  of  the  boundaries  and  to  the  number  of  material  types  that  make  up  the  object. 
Further,  differential  methods  can  model  extreme  changes  of  material  characteristics  with  ease  and  can 
handle  material  nonlinearity  efficiently.  The  performance  of  integral  methods,  on  the  other  hand, 
deteriorates  not  only  as  the  size  of  a  scatterer  grows  but  also  as  its  structure  becomes  more 
complicated.  Unlike  the  case  with  differential  methods,  the  size  of  integral  solution  matrices 
increases  with  the  addition  of  each  new  boundary  in  a  structure  of  given  size  and  modeling  nonlinear 
material  characteristics  with  integral  methods  is  difficult. 


In  the  final  analysis,  the  efficiency  of  any  method  for  solving  electromagnetic  scattering  problems 
must  be  judged  by  running  actual  computer  programs.  The  performance  of  any  algorithm  is  affected  by  the 
implementation  and  sometimes  synergisms  exist  in  implementations  that  are  difficult  to  quantify.  However, 
the  above  analysis  shows  that  differential  methods  are  clearly  favored  as  scattering  problems  become  larger 
and  more  complicated.  In  the  five  categories  examined  above,  where  there  is  a  difference,  differential  methods 
prove  to  be  superior  to  integral  methods.  In  particular,  in  terms  of  matrix  set  up  time,  matrix  solution  time, 
and  robustness  the  better  choice  is  the  differentia]  approach. 
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1.1.2  TIME  DOMAIN  VERSUS  FREQUENCY  DOMAIN 


A  second  issue  in  modeling  electromagnetic  scattering  from  dielectric  scatterers  is  the  choice  of 
analysis  in  the  time  or  the  frequency  domain.  With  continuous  wave  (CW)  signals,  the  logical  choice  is  to 
solve  for  the  electromagnetic  wave  in  the  frequency  domain  since  that  is  where  the  CW  signals  reside;  with 
pulsed  signals,  time  domain  methods  appear  to  be  the  most  appropriate  choice  since  the  signal  is  in  fact  a 
transient  in  time.  However,  as  discussed  more  fully  below,  time  domain  solutions  are  expensive  since  they 
require  many  time  steps  to  produce  acceptable  results.  For  large  problems  it  is  more  efficient  to  decompose 
the  transient  into  a  few  Fourier  components,  solve  the  problem  in  the  frequency  domain  at  these  few 
frequencies  and  then  reassemble  the  solution  in  the  time  domain.  In  addition  to  efficiency,  this  procedure  has 
the  advantage  that  the  dielectric  constant  of  many  materials  is  given  as  a  function  of  frequency;  with  time 
domain  solutions,  correct  material  property  values  are  difficult  to  determine  since  they  are  usually  given  as 
frequency  dependent  functions. 

Nevertheless,  all  published  work  for  modeling  three  dimensional  scattering  using  differential 
methods  is  in  the  time  domain.  Taflove  and  his  coworkers  [14,15]  have  employed  the  finite  difference  method 
in  the  time  domain,  while  Mei  and  his  coworkers  [12,13]  have  employed  the  finite  element  method  in  the 
time  domain.  Common  with  both  approaches  is  the  time  integration  formula  of  Yee  [16]  that  provides  a 
stable  time  integration  procedure  for  Maxwell's  equations.  The  Yee  integration  algorithm  is  the  key  to  time 
domain  solution  procedures:  it  is  stable  and  does  not  require  the  solution  of  matrix  equations.  As  we  will  see 
below,  stable  solution  procedures  for  three  dimensional  problems  in  the  frequency  domain  have  not  appeared 
in  the  literature  previously. 

The  major  advantages  of  the  finite  difference  time  domain  (FDTD)  method  is  that  it  is  extremely 
simple  and  that  it  is  unnecessary  to  store  matrix  entries.  Its  disadvantages  are  that  it  is  a  low  order  method  and 
that  it  employs  rectangular  meshes  with  minor  variations  in  element  size  and  shape.  With  the  finite  element 
method  in  the  frequency  domain  (FEFD),  although  a  matrix  equation  needs  to  be  solved,  the  accuracy  of  the 
solution  can  be  improved  by  using  higher  order  approximation  functions  and  by  employing  irregular  meshes. 

To  compare  the  relative  efficiency  of  the  FDTD  and  the  FEFD  methods,  we  will  make  an  order  of 
magnitude  comparison.  As  is  the  case  in  comparing  differential  and  integral  methods,  this  comparison  will 
neglect  all  factors  except  the  highest  order  trends  and  thus  will  serve  to  indicate  only  theoretical  limits  as  the 
problem  '  ize  grows.  To  begin,  the  number  elements  N  to  achieve  a  given  accuracy  in  three  dimensional 
problems  goes  as:  N  ■  h'^  where  h  is  the  mesh  size  used.  Since  the  FDTD  algorithm  is  based  on  linear  field 
interpolation,  the  order  of  accuracy  of  the  FDTD  algorithm  is  h^.  It  has  been  shown  [12-15]  that,  to  obtain 
accurate  solutions,  the  time  step  At  should  be  of  the  same  order  as  the  mesh  size:  At *•  h.  Thus,  the  number  of 
time  steps  to  obtain  the  solution  goes  as  h**.  The  the  total  complexity  to  obtain  a  FDTD  solution  goes  in 
the  limit  of  increasing  problem  size  as  N/h  ■»  N*-333. 

The  order  of  accuracy  of  the  finite  element  method  depends  on  the  order  of  polynomial  interpolation 
used.  As  a  reasonable  value  and  one  that  makes  the  calculation  relatively  easy,  we  will  use  cubic  order  finite 
elements  in  the  following.  The  order  of  accuracy  of  cubics  is  where  s  is  the  largest  edge  length  in  the 
mesh.  Assuming  a  uniform  mesh  (which  greatly  lowers  the  efficiency  of  finite  element  methods),  we  see  that. 


8 


to  achieve  the  same  accuracy,  the  FDTD  mesh  size  h  and  the  FEFD  mesh  size  s  are  related  by:  h  «  s^.  Thus 
the  number  of  finite  elements  M  required  to  achieve  the  same  accuracy  as  a  FDTD  mesh  of  N  elements  is  M  ■« 
VN  or  in  terms  of  mesh  size  M  «■  s"3  «  h'*A  In  the  FEFD  method,  one  must  solve  a  sparse  matrix  equation. 
As  discussed  before,  the  order  of  magnitude  work  required  to  do  this  increases  linearly  with  the  number  of 
elements  if  the  preconditioned  conjugate  gradient  algorithm  is  used  with  a  multigrid  mesh  and  goes  as  M1^  if 
a  single  level  mesh  is  employed.  The  net  result  is  that  the  total  complexity  if  the  FEFD  method  goes  either 
as  M «« N^-5  or  as  N^-6  depending  on  the  number  of  mesh  levels  employed.  In  either  case,  the  FEFD  method 
is  far  more  efficient  than  the  FDTD  method. 

It  must  be  emphasized  that  the  above  analysis  ignores  the  algorithmic  details  of  computing 
solutions  on  relatively  small  meshes  and  only  examines  global  trends  as  the  problem  size  grows.  Since  the 
FDTD  method  is  very  simple  while  the  FEFD  is  relatively  complicated,  it  may  be  that  for  small,  fixed-sized 
problems  the  FDTD  method  is  actually  faster.  The  above  analysis  shows,  however,  that,  as  the  problem  size 
gets  larger,  eventually  the  FEFD  method  will  be  the  most  efficient  Also,  the  FEFD  method  allows  the  use  of 
irregular  meshes  that  are  clearly  superior  for  modeling  problems  having  complicated  geometries. 


1.1.3  FINITE  DIFFERENCES  VERSUS  FINITE  ELEMENTS 


It  thus  appears  that  the  most  efficient  procedure  to  model  three  dimensional  electromagnetic 
scattering  is  to  use  differential  methods  in  the  frequency  domain.  There  remains,  however,  the  question  of 
which  differential  method  to  employ:  the  finite  difference  method  or  the  finite  element  method.  We  note  here 
that  the  differences  between  these  methods  is  often  just  semantics:  different  people  will  sometimes  call  the 
same  method  finite  elements  or  finite  differences  depending  on  their  viewpoint  For  example,  the  methods 
developed  by  Taflove  and  by  Mei  are  identical  except  that  Taflove  uses  a  rectangular  mesh  and  calls  it  finite 
differences  while  Mei  allows  the  mesh  to  be  distorted  and  calls  it  finite  elements.  Here  we  take  a  finite 
difference  method  to  be  one  that  is  derived  by  considering  the  difference  equations  formed  by  point  values  of 
the  differential  equation  without  reference  to  variational  calculus.  A  finite  element  method,  on  the  other  hand, 
is  taken  to  be  a  method  that  employs  basis  functions  to  approximate  the  field  in  a  variational  procedure  - 
either  a  Galerkin  method  or  the  Rayleigh-Ritz  method. 

From  the  point  of  view  of  accuracy  and  efficiency,  finite  element  methods  have  four  advantages  over 
finite  difference  methods.  The  first  of  these  advantages  is  the  fact  that  one  can  show  [17]  that  the  errors  in 
finite  element  methods  are  minimized  in  a  least  squares  sense  and  hence,  in  general,  will  be  more  accurate  than 
finite  difference  methods.  Second,  since  the  finite  elements  may  be  of  any  shape,  the  meshes  formed  in  finite 
element  methods  are  often  more  flexible  geometrically  than  the  meshes  formed  in  finite  difference  methods. 
Third,  it  is  relatively  easy  with  finite  element  methods  to  construct  elements  having  a  high  order  of 
interpolation;  high  order  formula  are  hard  to  construct  with  finite  difference  methods.  And  fourth,  finite 
element  methods  satisfy  certain  boundary  conditions  automatically  -  called  the  natural  boundary  conditions  - 
that  are  lacking  in  finite  difference  methods. 
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1.1.4  SCOPE  OF  THIS  PROJECT 


For  the  above  reasons,  the  direction  of  research  pursued  in  this  project  was  to  develop  finite  element 
methods  in  the  frequency  domain  for  the  solution  of  open  three  dimensional  electromagnetic  scattering 
problems.  This  was  an  ambitious  goal:  prior  to  the  work  in  this  project,  no  stable  finite  element  method  in 
the  frequency  domain  for  solving  even  closed  cavity  three  dimensional  electromagnetic  fields  was  known. 
Further,  methods  of  coupling  the  finite  element  solution  region  to  the  exterior  field  were  grossly  inefficient. 
Unfortunately,  due  to  budgetary  constraints  the  manpower  available  to  work  on  this  project  was  limited  and  it 
was  not  possible  to  develop  a  three  dimensional  scattering  program  during  the  course  of  this  project. 
However,  a  great  deal  of  useful  work  was  accomplished.  In  particular,  we  solved  the  two  bottleneck  issues  in 
applying  finite  element  analysis  to  this  area:  in  this  project  we  developed  and  tested  for  the  first  time  stable 
finite  element  methods  for  the  steady-state  three  dimensional  analysis  of  electromagnetic  fields,  and  we 
developed  the  transfinite  element  method  for  the  analysis  of  unbounded  problems  in  electromagnetics  and 
tested  this  method  in  two  dimensions.  Thus,  although  we  did  not  accomplish  the  original  project  goal  of 
solving  the  electromagnetic  scattering  from  large  three  dimensional  dielectric  objects,  we  are  proud  to  have 
developed  three  important  new  concepts  -  tangential  finite  elements,  C1  quadratic  finite  elements,  the 
transfinite  element  method,  as  well  as  the  lesser  but  useful  adaptive  spectral  response  modeling  procedure  - 
during  the  course  of  this  project.  It  should  be  noted  that  work  reported  here  has  been  invited  for  presentation  at 
two  conferences,  one  to  be  given  at  the  IEEE  APS/URSI  Symposium  to  be  held  in  Syracuse,  NY  in  June, 
1988  and  the  second  invited  paper  to  be  given  at  the  IEEE  Intermag  Conference  to  be  held  in  Vancouver,  BC 
in  July,  1988. 
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SECTION  2  FINITE  ELEMENT  FORMULATIONS 

2.1  THE  ZERO  DIVERGENCE  CONDITION 


Although  the  application  of  the  finite  element  method  to  three  dimensional  electromagnetic  field 
problems  appears  at  first  glance  to  be  straightforward,  this  in  fact  is  not  the  case.  As  described  more  fully  in 
Appendix  A,  the  application  of  the  standard  finite  element  method  to  three  dimensional  electromagnetics 
results  in  "spurious  modes”.  These  spurious  modes  make  the  solution  procedure  worthless.  It  should  be  noted 
that  the  problem  of  spurious  modes  is  not  limited  to  closed  cavity  problems;  finite  element  solutions  of  any 
electromagnetics  problem  -  open  or  closed  -  will  be  spurious  if  the  basic  approximation  procedure  is  incorrect. 

Although  we  now  know  that  it  is  not  true,  the  common  view  published  in  the  literature  [3-11]  is 
that  spurious  modes  are  a  result  of  not  enforcing  the  zero  divergence  of  the  electric  or  of  the  magnetic  field 
rigorously.  This  view  was  created  by  the  fact  that  spurious  modes  do  exhibit  a  non-zero  divergence.  To  solve 
for  the  electromagnetic  field,  the  finite  element  method  is  applied  to  the  vector  wave  equation.  While  the  zero 
divergence  of  the  field  is  a  consequence  of  this  equation,  zero  divergence  is  not  enforced  explicitly  in  the  usual 
finite  element  procedure.  Several  published  papers  [9-11]  attempt  to  eliminate  spurious  modes  by  attempting 
to  enforce  the  zero  divergence  condition  on  the  vector  field.  In  the  published  work,  this  is  achieved  only  in  a 
least-squares  sense  by  adding  a  penalty  term  proportional  to  the  square  of  the  divergence  to  the  variational 
principle 


At  the  start  of  this  project,  we  also  believed  that  the  way  to  eliminate  spurious  modes  was  to 
enforce  the  zero  divergence  condition  on  the  electromagnetic  field.  In  fact,  at  the  time  that  we  wrote  the 
proposal  for  this  project,  we  had  just  developed  a  numerical  procedure  by  which  the  zero  divergence  condition 
could  be  imposed  explicitly  in  the  finite  element  method.  This  procedure  is  published  in  reference  [18].  This 
procedure  did  exactly  what  we  believed  was  required  to  eliminate  spurious  modes:  it  enforced  the  zero 
divergence  of  the  field  exactly.  Consequently,  we  stated  in  our  proposal  that  we  would  employ  this  procedure 
to  generate  stable  solutions  for  three  dimensional  electromagnetics  problems;  we  began  work  on  this  project 
believing  that  this  approach  would  solve  the  instability  problem. 

Unfortunately,  however,  we  discovered  later  that  the  zero  divergence  condition  is  in  fact  unrelated  to 
the  problem  of  spurious  modes.  Enforcing  the  the  zero  divergence  condition  by  the  method  of  reference  [18] 
does  not  alter  the  fundamental  characteristic  of  standard  finite  element  methods:  the  solution  is  unstable  at 
certain  unpredictable  frequencies.  The  idea  that  we  had  worked  towards  -  namely  to  enforce  the  zero  divergence 
condition  exactly  on  the  vector  field  -  did  not  accomplish  the  desired  result.  The  failure  of  this  approach  to 
produce  stable  three  dimensional  field  solutions  was  both  disappointing  and  perplexing. 
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2.2  THE  ORIGIN  OF  SPURIOUS' MODES 


Fortunately,  a  publication  in  the  IEEE  Transactions  on  Microwave  Theory  and  Techniques  set  us  on 
the  path  the  would  lead  to  the  eventual  solution  of  the  spurious  mode  problem.  A  paper  by  Hano  [19]  showed 
that  spurious  modes  do  not  appear  in  two  dimensional  finite  element  solutions  if  rectangular  mixed 
constant-linear  elements  are  used  as  the  basis  functions  in  the  approximation.  While  rectangular  elements  are 
not  ideal  for  modeling  complicated  shapes,  Hano's  paper  did  allow  us  to  derive  the  key  principle  required  for 
stable  solutions  of  the  electromagnetic  field:  It  must  be  possible  to  express  the  basis  functions  used  in  the 
finite  element  method  as  the  gradient  of  a  scalar  for  this  method  to  provide  stable  solutions  of  the  vector  wave 
equation.  As  will  be  described  below,  we  have  used  this  principle  to  derive  several  different  methods  for 
computing  stable  solutions  to  electromagnetic  field  problems. 

To  understand  above  principle,  notice  that  the  curl  of  the  gradient  of  any  scalar  function  is  zero. 
Thus,  these  functions  are  nullvectors  of  the  curl  operator.  Since  the  operator  in  the  vector  wave  equation  is  the 
curl  operator  twice  over,  the  gradient  of  any  scalar  is  an  eigenvector  of  the  vector  wave  equation  with 
eigenvalue  zero.  Provided  that  the  finite  element  approximation  functions  can  be  expressed  as  the  gradient  of  a 
scalar,  the  approximation  functions  themselves  provide  a  basis  for  the  nullspace  of  the  curl  operator.  The 
eigenvalues  of  the  approximation  functions  will  thus  be  zero.  However,  if  the  approximation  functions  cannot 
be  expressed  in  the  required  form  then  their  eigenvalues  will  not  be  zero.  Hence,  the  eigenvalues  of  the 
nullspace,  which  should  all  be  zero,  will  be  shifted  to  non-zero  values.  This  is  the  source  of  spurious  modes: 
spurious  modes  are  approximations  to  the  nullvectors  of  the  vector  wave  equation  that  are  shifted  in  frequency 
so  much  that  they  interfere  with  the  real  physical  modes. 


2.3  STABLE  FINITE  ELEMENTS 


We  are  thus  able  to  explain  the  cause  and  the  origin  of  spurious  modes.  The  next  step  is  to  develop 
finite  elements  with  the  required  properties  to  avoid  them.  The  first  way  to  accomplish  this  is  to  extend 
Hano's  element  to  three  dimensions  and  to  higher  order  polynomials.  This  is  relatively  straightforward;  the 
mathematical  details  are  presented  in  Appendix  A.  As  shown  by  the  example  problems  in  this  Appendix, 
these  elements  work  exactly  as  predicted:  no  spurious  modes  occur  when  solving  three  dimensional  field 
problems  using  these  new  elements.  There  are  zero  eigenvalues  but  these  zero  eigenvalues  do  not  interfere 
with  the  physical  solutions.  Furthermore,  the  number  of  zero  eigenvalues  equals  exactly  the  dimension  of  the 
nullspace  of  the  curl  operator. 

While  the  success  of  the  rectangular  elements  in  Appendix  A  in  eliminating  spurious  modes  is 
encouraging,  it  is  not  sufficient.  Rectangular  elements  are  just  as  crude  for  approximating  complex  shapes  as 
are  finite  difference  methods.  Clearly,  triangular  finite  element  methods  that  satisfy  the  above  property  are 
required  if  we  are  to  achieve  the  efficiencies  hoped  for  at  the  start  of  this  project 
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Finding  triangular  finite  elements  that  represent  the  gradients  of  scalars  is  not  an  easy  task.  Suppose 
that  we  choose  the  scalar  to  be  an  ordinary  triangular  finite  element  complete  in  polynomials  of  order  n.  It  is 
easily  seen  that  the  derivative  of  a  complete  n'th  order  polynomial  is  a  complete  polynomial  of  order  (n  - 1). 
However,  since  (n  -  l)'st  order  finite  elements  are  basically  the  same  as  n'th  order  elements,  we  are  lead  to  the 
false  conclusion  that  the  components  of  the  gradient  can  be  expressed  in  terms  of  ordinary  triangular  finite 
elements.  We  know  that  this  conclusion  is  false:  ordinary  finite  elements  generate  spurious  modes  in  general! 

To  resolve  the  above  paradox,  we  need  to  consider  not  only  the  effect  of  the  gradient  on  the  inside  of 
the  element  but  also  its  effect  along  the  element  edges.  Ordinary  triangular  finite  elements  are  said  to  be 
continuous.  This  means  that  the  functions  in  these  approximation  are  continuous  but  that  their  derivatives  are 
not  continuous.  The  gradient  of  a  finite  element  will  therefore  be  a  vector  that  has  continuous  tangential 
components  but  discontinuous  normal  components.  Thus,  stable  triangular  finite  elements  are  obtained  by 
creating  elements  that  impose  continuity  of  the  tangential  component  of  the  vector  field  but  allow  the  normal 
component  to  be  discontinuous.  Such  finite  elements  are  indeed  possible;  Appendix  B  presents  details  of  the 
mathematical  derivation.  It  should  be  noted  that  although  it  is  possible  to  create  these  tangential  vector 
approximation  functions,  this  does  not  guarantee  that  they  will  provide  correct  solutions  to  the 
electromagnetic  field;  the  fact  that  they  do  is  a  consequence  of  the  natural  boundary  conditions  in  the 
variational  procedure.  The  full  details  of  the  synergism  between  tangential  elements  and  the  natural  boundary 
conditions  of  the  variational  procedure  are  given  in  Appendix  B. 

There  is,  however,  a  second  solution  to  the  above  paradox.  Let  us  repeat  the  above  discussion,  but 
replace  the  C®  scalar  with  a  Cl  scalar.  With  a  C1  scalar,  both  the  function  and  its  first  derivative  are 
continuous  so  that  both  the  tangential  and  normal  components  of  its  gradient  are  continuous.  The  components 
of  the  gradient  are  thus  ordinary  C°  finite  elements;  this  time,  however,  they  must  be  compatible  with  a  C1 
scalar  as  its  parent. 

The  lowest  order  polynomials  that  allow  C*  continuity  are  quadratics.  Surprisingly,  no  one  had 
developed  an  explicit  form  for  these  polynomials  prior  to  the  work  reported  here.  The  basic  idea  with  C* 
quadratics  is  that  it  is  impossible  to  generate  derivative  continuous  polynomials  one  element  at  a  time;  rather, 
one  must  subdivide  each  triangular  element  in  a  certain  way  into  six  subtriangles  in  order  to  achieve  C* 
continuity.  The  full  details  of  this  process  are  given  in  Appendix  C. 

Thus  the  C*  elements  developed  in  Appendix  C  provide  a  fascinating  irony.  It  turns  out  that 
ordinary  finite  elements  work  correctly  after  all,  but  only  if  they  are  placed  in  a  special  pattern.  If  an  arbitrary 
pattern  of  finite  elements  is  used  to  approximate  the  components  of  electric  or  magnetic  field  vector,  then  in 
general  the  nullspace  of  the  curl  operator  will  be  improperly  modeled  and  spurious  models  will  result.  If, 
however,  the  pattern  of  finite  elements  that  approximate  the  vector  field  is  such  that  a  C1  scalar  is  possible, 
then  no  spurious  modes  will  result.  Note  that  it  is  not  required  in  this  process  to  form  the  C*  scalar 
explicitly;  the  possibility  of  these  C*  scalars  is  sufficient  to  ensure  stability  of  the  electromagnetic  solution. 

Finally,  the  C*  elements  may  themselves  be  used  directly  to  generate  stable  finite  element 
solutions  of  electromagnetic  field  problems.  Since  the  electric  and  magnetic  fields  are  nondivergent,  one  may 
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define  a  vector  potential  that  provides  these  fields.  The  divergence  of  this  potential  function  is  arbitrary;  all 
that  is  required  is  that  its  curl  equal  the  field  variable.  Substituting  this  vector  potential  into  the  vector  wave 
equation  gives  a  new  equation  that  requires  continuity  of  the  potential  and  of  its  derivative.  The  C1  finite 
elements  of  Appendix  C  have  exactly  the  required  properties.  Consequently,  if  we  use  these  C1  elements  to 
approximate  each  component  of  the  vector  potential  function,  stable  finite  element  solutions  result.  Details  of 
this  procedure  are  presented  in  Appendix  D. 

Since  we  have  now  developed  not  one  but  three  different  stable  procedures  for  solving  three 
dimensional  electromagnetic  field  problems  using  triangular  finite  elements,  it  is  appropriate  to  ask:  Which  is 
the  best?  From  the  point  of  view  of  simplicity  of  formulation  and  ease  of  use,  the  answer  is  clear:  The 
tangential  finite  elements  of  Appendix  B  are  simpler  to  understand  and  to  program  than  the  C*  elements  of 
Appendix  C  or  the  potential  formulation  of  Appendix  D.  Fortunately,  this  simplicity  also  translates  into  a 
computational  advantage:  since  the  matrices  produced  by  tangential  elements  are  very  sparse,  while  the 
matrices  produced  by  the  Cl  elements  are  relatively  full,  tangential  elements  are  not  only  simpler  to  use  than 
the  Cl  elements,  they  are  faster. 


3.  THE  INFINITE  SOLUTION  REGION 


3.1  EXTERIOR  FIELD  MODELS 


A  basic  characteristic  of  scattering  problems  is  that  the  solution  region  is  infinite.  A  plane  wave 
originates  at  infinity  and  travels  through  free  s^ace,  impinging  on  an  object  that  scatters  the  wave  in  various 
directions  back  towards  infinity.  In  order  to  solve  this  problem  with  finite  elements,  the  problem  region  must 
be  divided  into  two  parts:  (1)  an  interior  region  surrounded  by  a  boundary  called  a  "picture  frame"  that 
completely  encloses  the  scatterer,  and  (2)  an  exterior  region  that  extends  from  the  picture  frame  to  infinity. 

The  problem  of  modeling  the  exterior  field  region  is  difficult  because  the  solution  region 
encompasses  an  infinite  number  of  wavelengths.  Furthermore,  while  the  energy  in  the  scattered  wave  is  finite, 
the  energy  in  the  incident  plane  wave  is  infinite;  this  fact  complicates  the  use  of  energy-based  solution 
methods.  Despite  these  difficulties,  however,  a  variety  of  methods  have  been  developed  to  solve  scattering 
problems  by  finite  element  methods.  Exterior  field  formulations  may  be  divided  into  three  groups:  those 
based  on  integral  methods,  those  base  on  differential  methods,  and  those  based  on  modal  analysis.  In  the 
following,  we  list  the  essential  characteristics  of  each  approach. 


3.1.1  BOUNDARY  INTEGRAL  AND  BOUNDARY  ELEMENT  METHODS 


The  oldest  way  to  model  open  boundary  problems  with  finite  element  methods  is  simply  to  resort 
to  integral  methods  to  handle  the  infinite  exterior  region.  Such  solution  procedures  are  usually  called  boundary 
integral  methods.  A  common  variant  of  boundary  integral  methods  is  the  boundary  element  method  in  which 
the  boundary  integral  is  discretized  by  using  the  piecewise  polynomials  first  employed  in  the  finite  element 
method  [20].  With  these  methods,  one  has  to  compute  the  values  of  an  integral  that  is  singular  on  the 
boundary.  A  way  to  circumvent  this  difficulty  is  to  employ  two  picture  frames,  one  for  the  equivalent 
boundary  integral  sources,  another  for  the  observation  points.  This  procedure  was  employed  early  in  the 
history  of  numerical  methods  by  Silvester  [21]  and  by  Wexler  [22].  The  separate  source-observation  approach 
has  been  converted  into  a  pure  integral  equation  form  by  Ludwig  recently  [23]. 

The  advantage  of  boundary  integral  and  boundary  element  methods  is  that  they  are  very  well  known. 
Their  disadvantage  is  that  they  are  relatively  expensive  because  it  is  necessary  to  compute  the  required 
integrals  of  the  Green's  functions  and  because  the  boundary  elements  produce  a  full  matrix  on  the  picture 
frame  boundary. 
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3.1.2  BALLOONING,  CLONING  AND  INFINITESIMAL  SCALING 


Ballooning  is  an  entirely  different  approach  to  modeling  the  exterior  field  region.  It  was  invented  by 
Silvester  [24]  in  1974  and  has  the  advantage  of  being  a  pure  differential  approach.  The  idea  is  the  following:  A 
single  layer  of  finite  elements  is  created  outside  the  picture  frame  boundary.  This  layer  has  elements  placed 
such  that  the  radial  edges  all  converge  to  a  single  point  called  the  star  point.  By  "blowing"  this  layer  of 
elements  up  so  that  the  inner  boundary  of  the  new  layer  coincides  with  the  outer  boundary  of  the  original 
layer,  a  second  layer  of  elements  is  formed  that  is  similar  in  a  geometrical  hence  to  the  original  layer.  These 
two  layers  are  then  joined  and  the  nodes  on  the  common  boundary  are  eliminated.  The  resulting  double  layer 
has  the  same  number  of  unknowns  as  the  original  single  layer  but  is  more  than  twice  as  large.  It  is  then  also 
blown  up,  the  two  double  layers  joined  and  the  in  between  common  nodes  eliminated.  This  process  is  repeated 
until  it  approximates  the  infinite  region  as  closely  as  desired.  (In  practice,  seven  or  eight  iterations  are  usually 
sufficient).  The  exterior  field  approximation  is  then  combined  with  the  interior  field  approximation  and  the 
problem  solved. 

The  original  approach  of  Silvester  could  only  be  used  to  model  Lapiacean  regions.  This  is  because 
in  two  and  three  dimensions  the  size  of  the  finite  elements  in  each  layer  grows  larger  and  larger.  In  a 
Lapiacean  region  this  is  desirable  because  the  solution  becomes  smoother  with  increasing  distance  from  the 
object,  but  this  is  inappropriate  with  wave  problems  where  the  distance  between  waves  is  invariant  with 
distance  from  the  scatterer.  The  ballooning  approach  was  extended  to  wave  propagation  by  Das  Gupta  and 
renamed  cloning  [25].  In  cloning,  only  a  single  layer  of  finite  elements  are  required,  but  the  property  of 
geometric  similarity  is  used  to  create  an  eigenvalue  equation  that  represents  the  entire  exterior  field  region. 

A  third  variant  of  ballooning  was  developed  by  Hurowitz  [26].  In  this  approach,  called  infinitesimal 
scaling,  the  single  layer  of  finite  elements  is  made  infinitely  thin.  In  this  limit,  the  exterior  field  problem  is 
converted  into  an  ordinary  differential  equation.  Unfortunately,  this  differential  equation  is  nonlinear  and 
expensive  to  solve. 


3.1.3  MODAL  ANALYSIS  METHODS 


The  third  procedure  for  solving  open  scattering  problems  is  semi-analytical.  In  this  case,  the 
electromagnetic  field  in  the  exterior  region  is  expressed  as  a  linear  combination  of  orthogonal  modes  and  the 
field  in  the  interior  region  is  expressed  as  a  linear  combination  of  finite  elements.  The  entire  solution  region 
is  solved  by  tying  these  two  approximations  together  with  continuity  conditions  along  the  picture  frame 
boundary. 


The  modal  approximation  method  for  electromagnetic  scattering  was  once  again  first  proposed  by 
Silvester  for  one  dimensional  propagation  to  infinity  [27].  Electromagnetic  scattering  in  one  dimension 
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assumes  that  the  scattered  field  is  confined  to  waveguides  that  are  connected  to  a  two  dimensional  waveguide 
junction.  In  Silvester's  approach,  the  electromagnetic  field  in  the  waveguides  are  expressed  in  terms  of  normal 
modes  and  the  normal  modes  in  the  waveguide  junction  are  solved  by  using  the  finite  element  method. 
Continuity  of  the  field  between  the  two  sets  of  modes  is  then  used  to  determine  the  total  field. 

The  modal  approach  to  two  dimensional  scattering  was  developed  by  Mei  [12,13]  and  called  the 
unimoment  method.  In  Mei's  approach,  the  electromagnetic  field  in  the  exterior  region  is  expressed  in  terms 
of  the  normal  modes  for  this  open  region.  Each  mode  is  then  used  to  set  up  a  finite  element  problem  that  is 
solved  for  the  interior  field.  The  total  field  is  given  by  a  linear  combination  of  these  finite  element  solutions. 
An  unfortunate  aspect  of  this  approach  is  that  continuity  must  be  imposed  on  both  the  field  and  its  derivative 
along  the  picture  frame  boundary. 

While  both  of  these  modal  solution  procedures  work,  they  are  both  expensive.  Silvester's  approach 
requires  that  the  normal  modes  of  the  finite  element  region  be  computed;  finding  these  eigenfunctions  is  a 
very  costly  process.  Mei's  approach  requires  that  n  finite  element  problems  be  solved  where  n  is  the  number 
of  modes  used  in  the  exterior  region.  Further,  Mei's  method  results  in  non-symmetric  matrices  with  which  the 
standard  preconditioned  conjugate  gradient  algorithm  doesn't  work. 

3.2  PRELIMINARY  INVESTIGATIONS 


At  the  start  of  this  project,  we  planned  to  use  the  cloning  algorithm  to  model  the  field  in  the 
exterior  region.  It  promised  to  be  fast  and  accurate,  and  was  simple  to  program.  As  stated  above,  in  cloning  a 
single  layer  of  finite  elements  is  created  on  the  surface  of  a  picture  frame  boundary.  The  entire  exterior  field 
description  is  extracted  from  this  single  layer  in  cloning,  while  in  ballooning  layer  after  layer  of  similar, 
larger  elements  are  added  to  the  first  The  idea  is  that  since  all  of  the  field  information  is  in  the  single  layer  - 
after  all  it  is  just  replicated  in  ballooning  to  form  the  other  layers  -  it  should  be  possible  to  obtain  the 
essential  characteristics  of  the  exterior  field  from  this  one  layer.  The  way  this  was  accomplished  by  Das  Gupta 
is  to  assume  that  the  field  falls  off  as  1/R  from  the  inside  to  the  outside  of  this  layer  imposing  this  condition 
on  the  single  layer  finite  element  model  results  in  an  eigenvalue  equation  for  the  field  in  the  single  layer  and 
hence  for  the  full  problem.  However,  in  using  the  algorithm  to  solve  some  test  problems,  we  were  unable  to 
get  the  correct  answers.  We  believe  now  that  decay  rate  of  1/R  is  correct  only  far  from  the  scatterer  and  does 
not  apply  to  a  layer  of  elements  that  are  close  to  the  scattering  object  To  make  the  procedure  valid,  we  must 
move  the  picture  frame  far  enough  from  the  scatterer  to  make  the  1/R  decay  rate  valid;  however,  at  that 
distance  there  is  no  need  to  compute  the  exterior  field  any  further  -  one  has  already  obtained  the  far  field 
pattern.  Consequently,  we  abandoned  the  idea  of  using  the  cloning  algorithm  of  Das  Gupta  to  compute  the 
exterior  field. 

The  next  approach  we  used  to  solve  the  exterior  field  problem  was  to  modify  the  ballooning 
algorithm  in  such  a  way  that  the  size  of  the  elements  in  the  replicated  layers  stays  constant  The  purpose  of 
this,  of  course,  is  to  ensure  that  if  the  first  layer  of  finite  elements  is  small  enough  to  approximate  the  wave 
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correctly,  then  all  other  layers  of  elements  are  small  enough  as  well.  To  accomplish  this  objective,  we  had  to 
abandon  the  idea  of  using  a  radial  star  point  to  generate  the  elements;  instead,  we  devised  various  rectangular 
patterns  of  elements  with  the  desired  property.  However,  it  became  apparent  in  working  with  these  patterns  in 
two  dimensions  that  this  approach  was  quite  unwieldy;  it  would  be  impossible  to  use  in  three  dimensions. 
Thus,  we  were  forced  to  abandon  the  ballooning  approach  as  well. 

Fortunately,  we  next  turned  to  the  modal  approach  to  model  the  exterior  field.  After  a  few  false 
starts,  we  discovered  an  excellent  way  to  model  the  exterior  field  in  scattering  problems.  We  call  the  new 
approach  the  transfinite  element  method,  and  it  is  described  next. 


3.3  THE  TRANSFINITE  ELEMENT  METHOD 


The  basic  idea  of  the  transfinite  element  method  is  to  use  the  modal  expansion  functions  that 
represent  the  exterior  field  as  approximation  functions  in  the  variation  procedure.  In  effect,  the  modes  form 
semi-infinite  finite  element  basis  functions  for  the  semi-infinite  region.  Full  details  of  the  transfinite  element 
method  along  with  some  examples  of  solutions  of  two  dimensional  scattering  problems  are  given  in 
Appendix  E. 

There  are  many  advantages  of  the  transfinite  element  method  compared  with  other  procedures: 

1.  It  is  a  direct  method.  Since  the  modal  functions  in  the  transfinite  element  method  are  used  in 
exactly  the  same  way  as  finite  element  functions,  this  method  gives  the  solution  of  scattering 
problems  directly  in  one  step.  In  the  unimoment  method,  modal  functions  are  used  to  match  the 
exterior  and  interior  field  solutions  but  are  not  used  directly  in  the  variational  procedure. 
Consequently,  many  finite  element  solutions  must  be  computed  to  obtain  a  single  field  solution 
with  the  unimoment  method,  while  only  one  finite  element  solution  is  required  for  each  frequency 
with  the  transfinite  element  method. 

2.  It  is  a  symmetric  method.  The  transfinite  element  method  produces  symmetric  matrix 
equations  that  are  easily  solved  by  the  preconditioned  conjugate  gradient  algorithm.  With  large 
problems,  this  is  orders  of  magnitude  faster  than  solving  non-symmetric  matrix  equations  by  other 
sparse  matrix  techniques. 

3.  It  is  a  simple  method.  Only  one  solution  is  required  and  only  continuity  of  the  field  is  required 
in  the  transfinite  element  method.  Methods  that  require  multiple  solutions  and  continuity  of  both 
the  field  and  its  derivative  are  more  difficult  to  program  and  are  less  accurate  as  well. 

4.  It  is  an  efficient  method.  With  almost  all  exterior  boundary  models,  the  number  of  variables 
corresponding  to  the  exterior  field  is  equal  to  the  number  of  nodes  on  the  picture  frame  boundary. 
This  is  very  costly:  there  are  many  picture  frame  boundary  nodes  and  all  of  these  nodes  are 
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interconnected  so  that  the  corresponding  matrix  elements  are  full.  In  the  transfinite  element  method, 
the  number  of  variables  corresponding  to  the  exterior  field  is  equal  to  the  number  of  modes  -  not 
nodes  -  used  to  approximate  the  exterior  field.  This  is  a  relatively  small  number  -  usually  on  the 
order  of  ten  or  twenty  -  compared  to  the  hundreds  or  even  thousands  of  nodes  on  the  boundary. 
Thus,  the  transfinite  element  method  provides  a  relatively  efficient  procedure  for  approximating  the 
exterior  field. 

In  every  respect,  the  transfinite  element  method  is  better  than  the  existing  alternatives.  By  making 
the  exterior  field  model  simple  and  efficient,  the  solution  of  large  scattering  problems  becomes  possible.  In 
fact,  the  transfinite  element  method  is  so  good  and  so  simple  that  one  is  tempted  to  ask:  Why  wasn't  the 
transfinite  element  method  invented  years  ago?  The  answer  seems  to  be  that  very  few  people  have  ever  worked 
on  differential  methods  to  model  scattering  problems,  and  the  few  that  have  worked  in  this  area  have  pursued 
other  directions  of  research. 
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4.  SPECTRAL  METHODS 

As  noted  in  the  Introduction,  frequency  domain  methods  are  more  efficient  for  continuous  wave 
modeling  than  are  time  domain  methods.  However,  many  problems  require  that  the  response  to 
electromagnetic  scattering  be  computed  over  a  wide  frequency  range.  If  done  by  brute  force,  this  still  may 
require  considerable  computer  time.  Often,  a  hundred  or  more  frequency  points  are  required  implying  that  a 
hundred  or  more  finite  element  problems  must  be  solved. 

In  an  effort  to  decrease  the  computational  requirements,  we  have  pursued  the  following  idea: 
Suppose  that  we  compute  the  field  problem  at  just  a  few  frequencies  and  then  use  these  frequencies  to  give  the 
solution  everywhere  within  the  band  of  interest.  This  linear  combination  of  solutions  will  provide  a  good 
approximation  to  the  field  at  other  frequencies  provided  that  the  solution  frequencies  account  for  all  of  the 
significant  events  in  the  scattering.  The  question  is  thus:  How  can  one  determine  the  proper  subset  of 
frequencies  in  a  frequency  range  that  gives  the  correct  overall  solution? 

To  answer  this  question  we  have  developed  an  adaptive  frequency  response  modeling  procedure. 
Essentially,  the  idea  is  a  simple  one:  First  solve  the  problem  twice  -  once  at  a  low  frequency  and  a  second 
time  at  a  high  frequency  -  and  then  evaluate  the  residual  of  the  linear  combination  solution  at  one  hundred 
points  in  between.  Then  solve  the  problem  again  at  the  frequency  that  gives  the  highest  residual,  form  the 
linear  combination  solution  with  all  three  solutions,  compute  the  residuals,  and  repeat  the  process  over  and 
over.  After  six  or  seven  iterations,  the  residuals  at  all  one  hundred  points  are  usually  within  acceptable  limits 
so  that  the  procedure  stops;  instead  of  computing  a  hundred  finite  element  solutions,  we  have  obtained  the 
same  accuracy  over  the  full  frequency  range  with  only  six  or  seven  finite  element  solution.  Details  of  the 
adaptive  spectral  response  modeling  procedure  are  given  in  Appendix  F.  Although  the  application  in  this 
Appendix  is  for  modeling  one  dimensional  electromagnetic  scattering,  no  change  is  required  to  use  this 
procedure  in  modeling  two  and  three  dimensional  scattering. 
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5.  CONCLUSIONS 

This  report  provides  several  breakthroughs  in  the  computation  of  electromagnetic  fields.  For  the 
first  time  stable  finite  element  solutions  of  three  dimensional  electromagnetics  problems  have  been  obtained. 
Further,  a  new  method  of  solving  two  dimensional  scattering  problems  called  the  transfinite  element  method 
has  been  developed  that  is  much  more  efficient  than  the  existing  alternatives. 

With  regard  to  the  issue  of  stability,  until  now,  no  one  could  be  certain  with  a  finite  element 
solution  whether  the  computed  solution  corresponded  to  a  real  physical  field  or  was  nonsense.  We  have, 
however,  explained  in  this  report  the  root  cause  of  these  spurious  solutions.  We  have  also  developed  three 
different  procedures  to  eliminate  them:  (1)  use  mixed-order  rectangular  elements,  (2)  use  edge-based  tangential 
continuity  elements,  or  (3)  use  derivative  continuous  C^  elements.  Although  the  mixed-order  rectangular 
elements  were  the  first  to  be  developed,  they  are  perhaps  the  least  useful  because  of  the  inflexibility  of 
rectangles  in  modeling  complex  shapes.  The  edge-based  and  the  C*  elements  are  both  tetrahedral  and  do  not 
suffer  from  this  geometric  limitation.  Of  these  two  elements,  the  edge-based  element  is  the  most  useful  since 
it  is  simpler  and  more  efficient  than  the  C1  element. 

With  regard  to  the  problem  of  efficiency,  two  major  developments  are  reported:  First,  for  the  first 
time  finite  element  solutions  of  open  electromagnetic  scattering  problems  have  been  obtained  by  using  a 
one-step  process  that  we  have  called  the  transfinite  element  method.  Prior  to  the  transfinite  element  method, 
scattering  problems  required  that  many  time  steps  or  many  modes  be  used  to  produce  each  finite  element 
solution.  The  one  step  transfinite  element  method  is  many  orders  of  magnitude  more  efficient  for  solving 
scattering  problems  than  the  existing  alternatives.  In  addition,  it  is  more  accurate  since  it  does  not  require  that 
derivatives  of  the  solution  be  computed  numerically  but  relies  on  variational  principles  to  set  derivative 
values. 


The  second  major  efficiency  improvement  reported  here  is  the  creation  of  an  adaptive  algorithm  to 
compute  the  scattered  fields  over  a  range  of  frequencies.  Heretofore,  the  computation  of  electromagnetic 
scattering  over  a  frequency  range  required  that  it  be  performed  perhaps  a  hundred  different  times;  we  have 
shown  that,  by  adaptively  selecting  the  solution  frequencies,  the  same  accuracy  throughout  the  frequency  range 
of  interest  can  be  obtained  by  using  only  six  or  seven  properly  chosen  frequencies.  Thus,  for  scattering 
problems  in  which  the  response  over  a  range  of  frequencies  is  desired,  the  adaptive  algorithm  provides  an  order 
of  magnitude  improvement  in  the  efficiency  of  the  solution  process.  The  algorithm  was  applied  in  this  report 
to  the  solution  of  the  one-dimensional  scattering  problems  encountered  in  waveguide  junctions;  extension  to 
two  and  three-dimensional  scatterers  is  straightforward. 

Unfortunately,  while  considerable  progress  has  been  reported,  there  was  insufficient  time  to  put  all 
of  the  pieces  together  and  develop  a  working  computer  code  to  solve  three-dimensional  scattering  problems. 
Although  we  pursued  this  research  in  a  logical,  methodical  manner,  one  graduate  student  working  for  thirty 
months  could  not  complete  die  desired  three  dimensional  scattering  computer  program.  In  fact,  the  nominal 
’one’  graduate  student  supported  by  this  project  turned  out  to  be  three  different  students  in  succession  since, 
due  to  various  circumstances,  two  students  left  this  project  and  had  to  be  replaced.  However,  several  major 
difficult  theoretical  problems  have  been  solved.  The  fact  that  we  did  not  succeed  in  solving  bistatic  scattering 
problems  from  large  three  dimensional  absorbing  bodies  (kies  not  mean  that  the  approach  taken  here  is  not  a 
good  one. 
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APPENDIX  A  -  MIXED-ORDER  RECTANGULAR  FINITE  ELEMENTS 
FOR  THE  SOLUTION  OF  THREE-DIMENSIONAL  ELECTROMAGNETIC 
FIELDS 


ABSTRACT 

A  new  method  for  modeling  electromagnetic  waves  by  the  finite  element  method  is  presented  in  this 
appendix.  In  this  formulation,  different  orders  of  polynomials  are  used  to  approximate  the  three  different 
components  of  either  the  electric  or  the  magnetic  '.eld  vectors.  We  show  that  this  procedure  eliminates 
the  problem  of  spurious  nodes  that  has  plagued  previous  three-dimensional  finite  element  solutions.  The 
method  is  applied  to  find  the  electromagnetic  fields  in  homogeneous  and  dielectric-loaded  cavities. 

INTRODUCTION 

The  finite  element  method  is  often  advanced  as  a  useful  numerical  procedure  for  modeling  high- 
frequency  electromagnetic  wave  phenomena.  Applied  at  an  early  date  to  solve  homogeneous  waveguide 
problems,1  the  method  has  proved  to  be  extremely  accurate  and  reliable  for  these  problems.  However, 
when  the  finite  element  method  was  applied  to  the  study  of  inhomogeneous  waveguides,  to  three- 
dimensional  resonant  cavities  and  to  scattering  problems,  difficulties  in  the  form  of  ■spurious*  modes 
were  encountered.  In  the  context  of  numerical  modeling,  a  spurious  mode  is  defined  to  be  a  non-pi  ysical 
solution  of  the  electromagnetic  field  equations  that  is  computed  simultaneously  with  the  correct  physical 
solutions. 

Since  the  presence  of  a  spurious  mode  in  a  numerical  solution  can  destroy  the  validity  of  the  solution, 
much  effort  has  been  directed  at  reducing  or  eliminating  the  unwanted  behavior.  The  first  approach, 
originally  suggested  by  Konrad,2  is  to  enforce  the  electromagnetic  field  boundary  conditions  exactly  on 
the  finite  element  approximation  space.  This  procedure  has  been  used  by  Mabaya,  Lagasse  and 
Vandenbulke3  in  an  E2  -  Hz  formulation  and  by  Davies,  Fernandez  and  Philippon4  and  by  Rahman  and 
Davies  5  in  the  3-component  H  formulation;  all  three  papers  report  only  limited  success  in  eliminating 

ft 

spurious  modes  by  this  technique.  Recently,  Koshiba,  Hayata,  and  Suzuki  have  shown  that  rigorous 
enforcement  of  boundary  conditions  does  indeed  eliminate  spurious  modes  above  the  "air  line"  (i.e.  the 
line  #/ko  =  1  in  a  0/ko  plot)  but  does  not  work  in  general  below  the  air  line. 


The  second  approach  to  eliminating  spurious  modes  is  to  modify  the  functional  in  the  variational  prin¬ 
ciple  used  to  approximate  the  fields.  Working  independently,  Winkler  and  Davies7  and  Hara,  Wada, 
Fukasawa  and  Kikuchi8  have  recognized  that  the  spurious  modes  do  not  satisfy  the  zero  divergence  con¬ 
dition  on  the  electric  or  magnetic  field.  Both  references  suggest  adding  a  penalty  term  proportional  to  the 
norm  of  the  divergence  of  the  field  to  the  functional.  Unfortunately,  this  procedure  does  not  eliminate 
the  spurious  modes  completely.  However,  as  demonstrated  in  reference  8,  the  spurious  modes  are  not 
stable  with  respect  to  the  amount  of  penalty,  and  can  be  distinguished  from  correct  solutions  by  plotting 
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the  finite  element  eigenvalue  spectra  with  respect  to  the  penalty  parameter.  Of  course,  this  procedure  is 
highly  inefficient  and  cumbersome:  each  new  field  problem  must  be  solved  repeatedly  and  all  of  the 
eigenvalues  plotted  in  order  to  identify  the  correct  solution. 


The  related  procedure  to  reducing  the  number  of  spurious  modes  was  proposed  recently  by  Konrad.9 
References  7  and  8  are  based  on  minimizing  a  functional  derived  from  the  vector  wave  equation  with  the 
addition  of  a  zero-divergence  penalty  term.  Konrad  suggests  using  the  vector  Helmholtz  equation  instead 
and  finds  that  "...  a  great  number,  though  not  all  of  the  spurious  solutions  are  indeed  eliminated."  This 
result  is  not  surprising  since  the  variational  expressions  derived  for  the  vector  Helmholtz  equation  and  for 
the  vector  wave  equation  with  the  addition  of  a  unit  penalty  term  are  identical  in  the  case  of 
homogeneous  media. 


The  third  approach  to  eliminating  spurious  modes  in  finite  element  solutions  is  to  restrict  the  finite 
element  approximation  functions  to  lie  in  a  reduced  vector  space.  In  this  view,  spurious  modes  are  the 
result  of  using  improper  functions  in  the  variational  procedure.  To  ensure  that  only  correct  solutions  are 
generated,  one  must  employ  only  admissible  functions  in  the  finite  element  approximation.  This  is  the 
approach  taken  in  this  report. 


The  use  of  a  restricted  function  space  to  eliminate  spurious  modes  in  finite  element  analysis  was  first 
suggested  by  Hano.10  Hano  showed  that,  in  two  dimensional  problems,  spurious  modes  are  completely 
eliminated  by  using  combination  constant-linear  finite  element  approximation  functions.  These  functions 
have  two  interesting  properties:  (l)  their  divergence  is  identically  zero,  and  (2)  they  are  discontinuous  at 
element  boundaries. 


In  this  paper,  we  derive  a  set  of  restricted  finite  element  basis  functions  for  the  solution  of  three 
dimensional  electromagnetic  field  problems  and  show  that  only  correct,  physical  solutions  are  obtained 
with  the  new  elements.  As  in  the  case  in  reference  10,  the  functions  reported  here  employ  different  orders 
of  polynomial  approximation  in  different  directions  in  each  element.  However,  in  addition  to  being  three 
dimensional,  the  approximation  functions  presented  here  are  high-order,  are  continuous  across  element 
boundaries  and  are  not  necessarily  non-divergent. 

FORMULATION  OF  THE  EM  FIELD  EQUATIONS 

Electromagnetic  wave  propagation  is  governed  by  the  vector  wave  equations. 
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V  X  —  v 


v  X  i  V 


E  -  «rVE 


X  H  -  ?rk0*B 


(1) 


(2) 


where  /if  ud  <r  ve  the  relative  permeability  and  relative  permittivity  of  the  material,  respectively,  aad 
k*  “  w5/i/0-  At  the  interface  between  two  dielectrics,  the  tangential  components  of  the  electric  aad 
magnetic  fields  must  be  continuous 

lm  x  (E(1)  -  E(2))  -  0  (S) 


lm  X  (H^  -  H(2))  —  0  (4) 

while  the  normal  components  are  discontinuous  as  follows 

1.  •  (<4 EW  -  CjE^2^)  -  0  (5) 


1.  •  (^hW  -  -  0  (6) 

In  these  equations,  superscripts  (l)  aad  (2)  refer  to  media  1  aad  media  2,  respectively,  and  1B  represents 
the  unit  normal  to  the  boundary. 

In  the  remainder  of  this  paper,  we  shall  use  the  electric  field  E  as  the  unknown.  Obviously,  a  similar 
treatment  holds  for  H.  In  terms  of  E,  equation  (4)  and  (5)  become 

1.  X  (—  V  X  E(I>  -  —  V  X  E(*>)  -  0  (7) 

*  "l 

tm  e  (VxE(1)  -  VxE(J))  -  0  (8) 

We  must  therefore  solve  equation  (l)  subject  to  the  interface  conditions  (S),  (5),  (7)  and  (8). 

Konrad1  has  shown  that  the  Euler  equation  of  the  functional 

F-\  J  |V  X  B|>  -  «,*,*»*>«  <») 

is  equation  (1)  and  that  the  corresponding  natural  boundary  equation  (7)  is  the  natural  boundary 
condition  for  this  functional.  If  one  wants  to  compute  the  electric  field  inside  a  cavity  containing 
variable-permeability  materials,  one  needs  to  minimise  (9)  in  an  appropriate  function  space. 

While  the  solution  process  described  above  has  been  widely  reported  in  the  literature*’4"7'1®"13,  as  noted 
in  the  introduction,  there  are  serious  problems.  To  eliminate  these  problems,  we  need  to  define  the 
approximation  functions  for  E  carefully. 
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A  BASIS  FOR  CURL 

To  solve  for  E  vis  equation  (9),  we  must  find  a.  legitimate  approximation  for  the  curl  operator.  An 
operator  has  a  domain,  a  nullspace  and  a  range;  it  is  not  enough  —  as  had  been  done  in  the  past  —  to 
approximate  only  the  domain  of  the  curl  operation. 

The  nullspace  of  an  operator  is  defined  to  be  the  set  of  functions  that  produce  sero  when  the  operator 
acts  on  it 

!^[A)  —  {*:  Ax  —  0  for  oil  x)  (10) 

It  is  well  known  that  the  nullspace  of  the  curl  operator  is  provided  by  the  gradient  operator 

tycurl)  -  V  d  (11) 

where  d  is  an  arbitrary  scalar. 

Let  us  approximate  d  by  finite  element  basis  functions  a(x,y,t)  over  a  rectangular  parallelepiped 

d(*,V,r)  —  o(,W)(*,y,*)£  (12) 

The  polynomial  a  (x,y,x)  is  m’lh  order  in  the  x-direction,  n’th  in  y,  and  p'th  in  s.  Using  the 
Kronecker  matrix  product,  three-dimensional  interpolation  polynomials  for  brick-shaped  elements  are 
given  by 

S *wW,r)  -  ff<m}(*)02r<"\p)®S^\r)  (IS) 

where  o^m^  is  a  one-dimensional  m’th  order  interpolation  polynomial  and  the  symbol  ®  denotes  the 
Kronecker  product  defined  in  Appendix  AA 

Since  the  derivative  of  an  n’th  order  polynomial  is  (n-l)’st  order,  it  follows  that  the  nullvectors  of  the 
curl  operator  must  have  the  form 

E,(z,V,t)  -  or 

(a,*— Is) 

-  S’  \  (M) 

EJ,z,y,t)  -  srtw-l)^ 

We  may  write  this  in  the  compact  form 

E  -  ng 


where 


(15) 
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1  — 


a 


0 

0 


E  - 


E„ 


E 

i_ 


0  0 


VxE  —  0CE 


a 


The  curl  of  E  is  evaluated  as 

~0 

•d/d  > 

e/e  7 

VxE  — 

did  s 

0 

-B/dx 

j9!B  7 

B/dx 

0 

Substituting  (IS)  into  (19)  gives 

E_ 


E_ 


E 
L  *  4 


(16) 


(17) 


(18) 


(18) 


where 


P  - 


0  grfw-U#-!) 

0 


0 

0 


(») 


c  . 


u-". 


where  the  matrices  D{  are  called  differentiation  matrices  and  are  defined  in  the  Appendix. 

We  note  that  the  factorisation  in  equation  (19)  is  not  possible  if  the  same  order  of  polynomial  is  used  to 
approximate  Ej  in  all  three  directions. 
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COMPUTING  THE  MATRIX  ELEMENTS 

Substituting  equation*  (15)  and  (10)  into  aquation  (0)  and  mini  mixing  with  respect  to  the  coefficients  E 
raaulu  in  the  matrix  aquation 

C^ICC  E  -  kf M  E  (21) 

where  K  and  M  are  the  matrices 


K  -  f  in  (22) 

dn  (23) 

To  evaluate  the  matrix  element*  in  K  and  M,  it  is  sufficient  to  evaluate  the  integral 


Substituting  equation  (13)  into  (25)  gives 


(24) 


G  -  (1  /L'L'L') J (cT(f) 0 q'(0 0 S^(X))r  (<T(f) 0 3"(fl 0 Sl(X))d^dX  (25) 

where  L  ,  L  and  L  are  the  dimensions  of  the  brick  element,  and  (,  {  and  X  are  homogeneous  coordinates 

*  w  * 

in  the  brick.  This  can  be  converted  into 

G  -  (26) 

where 


i<o  -  yvwfy) 


X(2)  « 


1 

6 


T<5)  «. 


1 

SO 


&%)ds 


(27) 

(28) 

(20) 


To  evaluate  the  finite  element  coefficient  matrix  for  one  dement,  one  therefore  needs  to  form  the 
metrics  K  and  M  using  equations  (27)  and  (28),  evaluate  C  by  using  the  values  in  the  Appendix,  and  pre- 
and  post-multiply  K  by  CT  and  C,  respectively.  The  contribution  from  each  dement  is  combined  with 
the  other  dements  in  the  grid  to  form  a  large,  sparse  matrix  eigenvalue  problem.  This  eigenvalue 
problem  is  solved  by  using  established  techniques  for  the  eigenvalues  kQ*  and  eigenvectors  g. 
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COMPUTATIONAL  RESULTS 

A  computer  program  has  beta  developed  baaed  on  the  above  formulation  to  aolve  for  the 
electromagnetic  fields  in  resonant  cavities.  The  program  allows  finite  elements  of  mixed  orders  to  be 
assembled  and  solved  for  complex  three-dimensional  geometries. 

EMPTY  CUBIC  BOX 

The  solution  provided  by  the  present  method  is  first  tested  with  the  analytical  solution  of  an  empty 
cubic  cavity  of  1  m*.  The  cubic  box  is  divided  into  2x2x2  elements.  The  results  obtained  by  using 
constant-bilinear  (m  —  n  »  p  1),  linear-biquadratic  (m  «  n  «  p  2)  and  quadratic-bicubic  (m  •  n  « 
p  3)  elements  are  shown  in  Tables  1-3.  A  one-to-one  correspondence  exists  between  the  approximate 
eigenvalues  and  the  exact  ones.  Since  no  spurious  modes  are  produced,  the  approximate  eigenvalues  are 
easily  identified.  The  results  also  show  that  the  dominant  eigenvalue  is  approximated  reasonably  well  even 
with  a  small  number  of  elements. 

Solutions  obtained  by  using  high-order  elements  are  aeen  to  be  much  more  precise  than  is  the  case  with 
constant-bilinear  element  aolutions.  The  reason  for  the  increased  accuracy  is  twofold:  (l)  Higher-order 
polynomials  are  more  accurate  than  low-order  ones,  and  (2)  higher-order  elements  are  continuous  across 
element  boundaries  while  constant-bilinear  elements  are  not. 

DIELECTRIC-LOADED  BOX 

Our  computational  procedure  has  also  been  applied  to  aolve  the  dielectric-loaded  cavities  shown  in 
Figure  1  and  the  results  compared  with  those  obtained  by  other  methods.  In  this  case,  the  magnetic  field 
is  solved  subject  to  natural  boundary  conditions  1bX(VxH)  «■»  0  on  perfectly  conducting  walls  and 
Dirichlet  conditions  lnXH  ■  Ood  symmetry  planes.  Comparison  of  our  results  with  previous  solutions 
is  shown  in  Table  4,  where  the  cavities  are  divided  into  2x2x2  linear-biquadratic  elements.  Again,  the 
dominant  eigenvalue  is  approximated  well,  despite  the  fact  that  only  8  elements  are  employed. 

NUMBER  OF  ZERO  EIGENVALUES 

The  number  of  the  unknowns  and  the  number  of  sero  eigenvalues  computed  depends  not  only  on  the 
number  of  dements  employed  but  also  on  the  type  of  boundary  conditions  imposed.  There  are  two  kinds 
of  vector  Dirichlet  boundary  conditions:  One  is  a  tangential  boundary  where  the  tangential  components 
of  the  field  vector  are  act  to  sero;  the  other  other  is  normal  boundary  where  the  normal  component  of  the 
field  vector  is  set  to  sero.  Empirically,  we  find  that  the  number  of  sero  eigenvalues  for  a  box  divided  into 
a  x  b  x  c  dements  is  given  by 

number  sero  —  (tea  +  s  —  fiaX<*  b  +  s  —  nb)(t»c  +  a  —  ne)  —  1  (30) 

+  number  of  disjointed  tangential  boundaries 

where  t  ■  l  and  i  ™  1  for  constant  bilinear  dements,  t  ■■  1  and  s  «  2  for  linear  biquadratic  dements,  t 
“  2  and  s  —  2  for  quadratic-bicubic  dements,  and  a  a,  nb  and  nc  are  numbers  of  Dirichlet  boundary 
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conditions  imposed  on  the  walls  normal  to  the  x,  y  and  s  directions,  respectively.  Enforcing  boundary 
conditions  exactly  does  reduce  the  sise  of  global  matrix  and,  consequently,  the  number  of  sero  eigenvalues. 
We  have  not  yet  been  able  to  explain  the  empirical  formula  in  equation  (30)  by  theoretical  means. 

CONCLUSIONS 

Three-dimensional  electromagnetic  field  problems  may  be  solved  by  using  mixed  order  finite  elements  to 
approximate  the  vector  field.  In  this  approach,  a  consistent  numerical  approximation  is  made  for  the 
domain,  range,  and  null  spaces  of  the  curl  operator.  Such  a  consistent  approximation  in  the  solution  of 
the  vector  wave  equation  eliminates  the  problem  of  spurious  modes  that  had  plagued  previous  solution 
procedures. 


A  limitation  of  the  present  formulation  is  the  restriction  of  the  element  shapes  to  be  rectangular 
parallelepipeds.  For  problems  involving  complicated  shapes,  tetrahedral  or  isoparametric  elements  would 
be  preferred.  The  use  of  such  elements  with  mixed-order  polynomial  basis  functions  is  presently  under 
investigation. 
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APPENDIX  AA  -  RECTANGULAR  BASIS  FUNCTIONS 


To  tee  bow  approximation  functions  nre  formed  for  bricks,  first  consider  the  two-dimensional  ease.  An 
arbitrary  rectangular  reference  element  is  mapped  onto  the  unit  square  by  the  transformation 


(*  -  *rVl. 

(Al) 

(v  -  vr)/l. 

(A2) 

Two-dimensional  basis  functions  may  be  obtained  from  one-dimensional  functions  by  using  the 
Kronecker  matrix  product. 


Definition: 

If  A  is  an  nxm  matrix  and  B  is  a  pxq  matrix,  then  the  Kronecker  matrix  product  of  A  and  B  is  denoted 
by  A®  B  and  is  the  npxmq  matrix 


ailB 


(AS) 


Kronecker  matrix  products  satisfy  the  following  useful  identities14 
(0  A®(B+C)-A®B+A®C 

(it)  oA  ®  0B—a0{A  ®  B) 

(««)  AB&CD-4A&CKB&D) 

(to)  (A  ®  B)_1— A"1  ®  B~l 
(v)  (A®B)T-AT®Br 


(A4) 


In  general,  one-dimensional  n’th  order  interpolation  polynomials  are  defined  as 
.  .  “I1  (f-f.) 

~  5 

where  the  ^  are  the  interpolation  nodes.  Using  the  Kronecker  matrix  product,  two-dimensional 
interpolation  polynomials  are  given  by 

s*,’",)(f.e)  -  a**\f)®aW(0  (A5) 


Now  consider  approximating  an  arbitrary  function  4(f,()  in  terms  of  the  finite  dement  approximation 
functions.  We  may  write  this  as 

*f.O  -  Z{mX()*a}*Xf)  (AS) 
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where  ♦  is  an  m  by  n  matrix  of  values  at  the  interpolation  nodes.  Equation  (A7)  is  converted  into  the 
standard  matrix  form 

*?,£)  -  a(r,0£  m 

by  defining  a  vector  operation  called  uee  in  the  following  manner: 


VecA 


(AS) 


where  A.  is  the  i’th  column  of  the  matrix  A.  The  operator  vec  and  the  Kronecker  product  are  related  by 
the  identity 

vec  ABC  —  {(?  ®  A)  vec B  (AS) 

The  vec  of  a  scalar  is  simply  itself;  therefore,  taking  vec  of  both  sides  of  equation  (AS)  yields 


-  &HXs)  03^(f))  £ 

where 

£  —  vec* 


(A10) 

(All) 


Three-dimensional  finite  elements  are  generated  by  an  analogous  procedure  to  that  used  in  two- 
dimensions.  Approximation  functions  for  brick-shaped  elements  are  given  by  the  equation 
S<W)  _  S<"\f)®ff(",\e)®o^)(X) 

(A12) 


where  (f,£,X)  are  homogeneous  coordinates  in  the  brick. 
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DIFFERENTIAL  MATRICES 

One-dimensional  differentiation  nitnctt  in  defined  by  the  equation 

d{ 


(Bl) 


where  is  the  n  by  n+1  differentiation  matrix.  The  polynomials  3^*  1\f)  in  equation  (Bl)  are  of 
one  order  leas  than  that  of  o^f)  because  the  derivative  of  an  n’th  order  polynomial  is  (n-l)’st  order. 
Evaluating  both  sides  of  equation  (Bl)  at  the  (n-l)’st  order  interpolation  nodes  —  1,.  .« 

provides  the  elements  of  the  differentiation  matrix  as 

~7Tk  -  .<-■> 


£}.(*)  ■«£>.(*)  mm 


•J 


•J 


Performing  the  indicated  operations  with  the  equispaced  node  interpolation  polynomials  gives  the 


numerical  values 

d<1>  - 

M. 

1] 

dW  - 

•3 

4. 

f 

*1. 

i. 

•4. 

3. 

— 

05.5 

9. 

-4.5 

1. 

.125 

-3.575 

3.375 

-.125 

-1. 

4.5 

-9. 

-5.5 

(B3) 


Note  that  the  matrix  D  has  the  following  anti-symmetry  property 

D  (•)  —  _n(*) 

Dii  * 


(B4) 


To  extend  the  above  result  to  two-dimensions,  we  need  to  evaluate 

I  -  i  sr«-lw®s<-\o 

Of  Of 

-  a?<-1J(f)D(w,®s<,l\e) 

-  («<m- «s>  sr(»J(OKi>(m)  ®  /) 


(B5) 


Thus  we  find  that 
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d 


-.0  - 


(Be) 


where 

£>f  —  D<m)®J  (B7) 

Similar  results  bold  for  derivatives  in  the  (-direction  and  for  three-dimensional  elements. 
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Table  1.  Eigenvalues  of  a  homogeneous  cubic  cavity  obtained  by  using  constant-bilinear  elements, 
approx-eigenvalue  modes  erect-eigenvalue 
0  8.881784  e-16 

1  2.400000e+0l  (1,1, P)  19  7892 

2  2.400000e+0l  (1,0,1)  19.7892 

*  2.400000e+0l  (0,1,1)  19.7892 

4  3.600000e+0l  (1,1,1)  29.6088 

5  3.600000e+0l  (1,1,1)  29.6088 


Table  2.  Eigenvalues  obtained  by  using  linear-biquadratic  elements. 


approx-eigenvalue  modes  exact-eigenvalue 

0  .1.1146640-18 

I  -9.870122O-14 

5  -7.9714010-14 
8  -6.372680O-14 
4  -4.7739590-14 
8  -3.1752380-14 

6  -1.5765170-14 

7  2.220446o-16 

•  1 .988769o+01  (1,1,0) 

0  1 .988769o+01  (1,0,1) 

10  1.9887690+61  (0,1,1) 

II  2.9615470+61  (1,1,1) 

12  2.9615470+61  (1,1,1) 

13  4.9943850+61  (2,1,0) 

14  4.9943850+61  (2,0,1) 

15  4.9943850+61  (1,2,0) 

16  4.9943850+61  (1,0,2) 

17  4.9943850+61  (04,1) 

18  4.9943850+61  (0,14) 

19  5.9266240+61  (2,1,1) 

20  5.9266240+61  (1,2,1) 

21  5.9286240+61  (1,14) 

22  5.947009O+61  (2,1,1) 

23  5.947009e+6l  (14,1) 

24  5.947009O+61  (1,14) 

25  8  000000e+0l  (2,2,0) 

25  8.000000o+61  (2,0,2) 

26  t.OOOOOOo+Ol  (044) 


19.7392 
19.7392 
19.7392 
29.6088 
29.6068 
494480 
494480 
494480 
494480 
49.3480 
494480 
59.2176 
594176 
594176 
594176 
59.2176 
59.2176 
7$.«?8 
75  9J48 
78.9568 
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Table  3.  Eigenvalues  obtained  by  using  quadratic-bicubic  elements.  The 
number  of  zero  eigenvalues  is  equal  to  64. 


•pprax-tiftaralue 

•odm 

■K*et«ceiiTalu« 

•4  1.074181*401 

U.1.0) 

10.7002 

65  1.074181*401 

(1,0,1) 

10.7082 

06  1.074101*401 

(0.1.1) 

10.7082 

67  0.060881*401 

(1.1.1) 

08.6088 

68  0.060801*401 

(1,1,1) 

00.6088 

60  4.087085*401 

(2,1,0) 

40.0480 

70  4.087085*401 

(2,0,1) 

40.0480 

71  4.087085*401 

(1.2,0) 

48.0480 

70  4.087085*401 

(1,0,2) 

40.0480 

70  4.087085*401 

(0,2,1) 

48.0480 

74  4.887085*401 

(0,1,2) 

40.0480 

75  5.070327*401 

(2.1.1) 

68.2176 

76  6.078027*401 

(1.2,1) 

68.2176 

77  5.070327*401 

(1,1.2) 

58.0178 

78  6.870610*401 

(2,1,1) 

50.2176 

70  6.870610*401 

(1X1) 

68.2176 

60  6.073610*401 

(1  A*) 

68.2176 

61  8.000000*401 

(2,2,0) 

78.0568 

62  0.000000*401 

(*.0,2) 

78.0568 

60  8.000000*401 

(0.2,2) 

78.0568 
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Table  4.  Dominant  resonant  frequency  KQ  a  in  Figure  1. 
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APPENDIX  B  -  EDGE-BASED  VECTOR  FINITE  ELEMENTS  FOR 
SOLVING  THE  ELECTROMAGNETIC  WAVE  EQUATIONS 

ABSTRACT 

This  appendix  examines  the  requirements  on  finite  element  approximation  functions  that  must  be 
satisfied  to  provide  consistent  solutions  of  electromagnetic  wave  problems.  We  show  that  continuity  of 
the  tangential  component  of  the  electric  or  of  the  magnetic  field  is  required  in  the  variational  procedure 
but  that  the  normal  component  of  the  field  need  not  be  continuous.  We  present  a  new  triangular  finite 
element  that  provides  for  only  continuity  of  the  tangential  component  of  the  field.  In  this  element,  the 
tangential  components  of  the  field  represent  the  unknowns  to  be  determined  and  are  evaluated  separately 
on  each  side  of  the  element.  The  finite  element  analysis  based  on  this  new  element  leads  to  the  complete 
elimination  of  the  spurious  modes.  Wave  propagation  in  dielectric-loaded  waveguides  is  examined  to  test 
the  validity  of  the  new  vector  element.  Extension  to  three-dimensionals  is  made  by  requiring  that  the 
tangential  components  of  the  field  interpolate  at  the  edges  of  tetrahedra.  This  3-D  vector  element  is 
shown  to  model  three-dimensional  cavity  problems  correctly  without  the  presence  of  spurious  modes. 

INTRODUCTION 

While  the  finite  element  method  is  known  to  be  a  powerful  tool  in  the  analysis  of  many  physical  sys¬ 
tems,  its  application  in  electromagnetics  has  been  limited  by  the  presence  of  spurious  or  non-physical  solu¬ 
tions  in  the  computation  of  electromagnetic  wave  problems.  Several  investigators  have  attempted  to 
eliminate  these  spurious  nodes  by  either  employing  a  penalty  method  or  by  strongly  imposing  boundary 
conditions.  However,  all  of  these  attempts  have  only  reduced  the  appearance  of  the  spurious  modes  to  a 
lesser  degree.  Recently,  a  rectangular  mixed-order  element  has  been  developed  based  on  the  consistency  of 
the  orders  of  the  polynomials  used  among  the  three  components  of  the  field  vector  in  the  three  inde¬ 
pendent  directions  1,z.  With  these  elements,  the  problem  of  spurious  modes  is  eliminated.  However,  the 
elements  reported  are  rectangles  in  2-D  or  rectangular  parallelpipeds  in  3-D  and  hence  provide  only  a 
limited  degree  of  geometrical  flexibility.  An  extension  of  this  concept  to  triangular  or  tetrahedral  finite 
elements  is  required  to  model  many  problems  efficiently. 

The  approach  reported  in  reference  [2]  can  be  traced  to  the  mathematicians  Raviart  and  Thomas3  and 
to  Nedelec4  who  were  the  first  to  seek  vector  solutions  in  the  range  space  of  the  operator.  These  resear¬ 
chers  have  presented  a  vector  finite  element  that  approximates  either  the  normal  or  the  tangential  com¬ 
ponents  of  a  vector  field  and  have  established  also  the  convergence  of  the  new  elements.  However,  the 
elements  reported  in  [3]  and  [4]  are  not  interpolator  and  hence  are  difficult  to  apply. 


41 


In  this  paper,  we  take  a  different  point  of  view:  We  examine  the  variational  principle  corresponding 
to  the  Maxwell  equations  and  show  that  electromagnetic  field  problems  are  best  solved  not  by  using  three 
separate  scalar  approximations  for  the  three  components  of  the  field,  but  by  using  new  vector  elements 
that  require  only  the  continuity  of  the  tangential  component  of  the  field.  With  this  new  element,  the 
normal  component  of  the  field  is  allowed  to  jump  across  the  material  interfaces.  Further,  we  improve 
Nedelec’s  tangential  elements  by  writing  the  field  vector  explicitly  in  terms  of  point  values.  This  inter¬ 
polation  property  makes  the  process  of  imposing  the  necessary  continuity  conditions  between  vector  ele¬ 
ments  easy  to  program  and  provides  a  practical  and  extremely  reliable  approach  for  solving  electromag¬ 
netic  field  problems. 

VARIATIONAL  FORMULATION 
ELEMENT  MATRICES  S  AND  T 
THE  ASSEMBLY  OF  GLOBAL  MATRICES 

Variational  procedures  for  solving  electromagnetic  field  problems  by  the  finite  element  method  were 
originally  proposed  by  Konrad12  and  used  subsequently  by  others13-15.  While  the  essence  of  Konrad’s 
derivation  is  correct,  the  interface  between  elements  were  not  examined  in  detail.  This  oversight  has  led 
to  the  appearance  of  spurious  modes  in  finite  element  solutions.  In  the  following  we  examine  these  equa¬ 
tions  in  detail,  paying  particularly  close  attention  to  the  interface  requirements  on  the  approximation  sub¬ 
space. 


The  basic  equations  that  govern  electromagnetic  fields  in  source-free,  time-harmonic  regions  are  the 
Maxwell  equations 


V  X  £  =  -  jwfiH 


(1) 


v  x  7?  =  jwiE  . 


(2) 


We  note  that  since  the  divergence  of  the  curl  of  anything  is  zero,  Equations  (l)  and  (2)  imply  that 
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fiTt  =  o 


v  •  tE  —  o  . 


(3) 


(4) 


The  substitution  of  E  from  Equation  (2)  into  (1)  and  of  7?  from  (l)  into  (2)  yields  the  vector  wave 
equations 

l2 


1  „  Jr  _ 

Vx-VxE  =  —  7/ 


v  x-v  x7?=  —  7? 

c  f 


(5) 


(6) 


n  n 

where  k  =  ur/i€.  In  view  of  the  similarity  between  Equations  (5)  and  (6),  the  following  analysis  will  be 
presented  in  terms  of  the  vector  E;  the  parallel  development  for  the  vector  7?  is  obtained  by  making  the 
substitutions  E  — »  7?  and  n  — *  (. 

INTERFACE  CONDITIONS 

The  interface  conditions  that  must  be  satisfied  by  the  electric  field  E  are  that  its  tangential  com¬ 
ponent  is  continuous 


ln  ><  =  “»  ><  E2 


(7) 


and  that  the  normal  component  of  the  flux  E  =  tE  is  continuous 


*  <1^1  =  Xn  *  <2^2  ' 


(8) 


For  the  magnetic  field  77,  the  corresponding  interface  conditions  are 


*«  X  *1  =  ln  ><  ^2 


1  n  *  "l77  =  1  n  # 


(9) 

(10) 


Equation  (10)  may  be  expressed  in  terms  of  the  electric  field  by  using  Equation  (1).  The  result  is 
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r  .  v  X  E,  -  1,  .  v  X  E2  .  (11) 

Adopting  an  (n,r,X)  coordinate  system,  where  1  is  normal  to  the  interface  and  1  f  and  1  x  are  parallel  to 
it,  Equation  (l  1 )  becomes 


Notice  that  Equation  (12)  only  involves  the  tangential  components  and  the  tangential  derivates  of  E  . 
We  therefore  arrive  at  the  important  result: 

Continuity  of  the  normal  component  of  magnetic  flux  ~B  =  pH  is  ensured  by  setting  the  tan¬ 
gential  component  of  the  electric  field  E  to  be  continuous. 

There  is  of  course  the  corollary: 

Continuity  of  the  normal  component  of  electric  flux  25  =  t~E  is  ensured  by  setting  the  tan¬ 
gential  component  of  the  magnetic  field  7/  to  be  continuous. 

It  follows  from  the  above  that  if  we  impose  the  continuity  of  the  tangential  components  of  E  and  7}  , 
Equations  (7)  and  (9),  directly  then  continuity  of  the  normal  components  of  the  fluxes,  Equations  (8)  and 
(10),  are  automatically  ensured. 

In  an  2?-field  solution  procedure,  we  must  express  the  tangential  component  of  the  magnetic  field  in 
terms  of  the  electric  field.  Using  Equation  (l),  the  result  is 

r»  X  7  (V  X  -  "n  X  7  (V  x  E2)  •  <13> 

h2 

Evaluating  this  by  components  yields  the  two  equations 


(15) 


Provided  that  Equations  (14)  and  (15)  are  satisfied,  Equation  (8)  will  be  satisfied. 
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THE  ENERGY  FUNCTIONAL 


Consider  the  functional 


f(E)  =  [  i  (v  x  E?<m  -  f  —  tfdn  . 

J  n  J  P 


(16) 


We  shall  show  that  this  functional  may  be  used  to  solve  for  the  electric  field  in  wave  problems  provided 
that  the  tangential  component  ojE  is  made  continuous. 

To  minimize  f\E),  let  E  =  E 4-  f  £  where  E a  is  the  exact  solution  of  the  wave  Equation  (6),  e  is 
a  number,  and  £  is  an  arbitrary  vector  function.  The  first  variation  of  F[E)  is 

af(E„  +  <{) 

,f,E)  -  - - -  Uo 


-  /  i  (V  X  a  .  (V  X  £ Jin  -  f  J  (  .  -Ejn 


(17) 


By  integrating  the  vector  identity 


V  •  (a  xF)  =  (VXa)*6  —  a  •  (V  X  F) 


(18) 


we  find  that 


£  (a  X  F)  •  dS  =  J  |(V  X  a )  •  6  —  a  •  (V  X  6  ) 


<m 


(19) 


Thus  (17)  is  converted  to 

6F\A)  =  /  T  • 


1  __  kl  _ 

v  x  -  V  x  E - E 

(t  a  (i  e 


+ 


/  X  J  v  X  •  Ts- 


(20) 


Since  Eg  satisfies  Equation  (6)  exactly,  the  volume  integral  in  (20)  drops  out,  leaving 


6F\E)  =  J>  ^  X  —  V  X  =  —juj){£  X  7?) 


7s 


(21) 
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Finally,  setting  the  first  variation  of  F(E)  equal  to  zero  yields 

f  (I  X  7?)  •  7s=  o  . 

Equation  (22)  provides  the  natural  boundary  conditions  for  the  functional. 
NATURAL  BOUNDARY  CONDITIONS 
Evaluating  Equation  (22)  by  components  yields 

/  -  0 

In  terms  of  the  electric  field  E  this  is 


1 

c 

A. 

dE 

n\ 

I 

c 

dE 
(  n 

s 

V 

^  dn 

dr) 

x  V 

^  d\ 

dn  J  ' 

On  exterior  boundaries,  if  we  set  the  tangential  components  of  E  equal  to  a  given  value,  then  £  and  £x 
are  zero  and  Equation  (24)  is  satisfied.  On  the  other  hand,  along  boundaries  where  we  leave  (  and  £x  to 

be  arbitrary,  we  will  automatically  obtain  a  solution  for  E  that  satisfies 


X  H  _  0 


On  interior  boundaries,  if  we  impose  continuity  of  the  tangential  component  of  E  then  we  will  have  that 

<J,  -  tJk 

Mi  “  Ma  ■  (« 


For  boundaries,  Equation  (24)  yields 

/•  , ,  JE.  dE. 


dE  dE_ 


[  [_L  (—1  _  _2\  _  L  fill  _  n\ 

J  nl  (an  dr  Ji  jij  (dn  dr  J 


dE  dE, 


dE  dE ' 


/I  ULJ.  N  1  UCj  . 

k  ( «r  -  srl  -  Z  («r  ‘  -  0 
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Since  and  are  arbitrary  on  interior  boundaries,  Equation  (27)  yields  exactly  the  correct  interface 
conditions  (14)  and  (15). 

We  note  that  the  continuity  of  the  normal  component  of  E  is  not  required  in  the  above  formulation. 

In  fact,  imposing  continuity  on  the  normal  component  of  E  is  wrong  because  1  •  E  is  discontinuous  at 

material  boundaries.  The  beauty  of  the  procedure  presented  in  this  report  is  derived  from  this  fact:  if 
one  allows  the  normal  component  of  the  electric  field  to  be  discontinuous  imposing  continuity  of  only 
tangential  component,  then  the  amount  of  discontinuity  in  the  normal  component  is  automatically  given 
by  the  natural  boundary  conditions  in  the  variational  principle. 

TANGENTIAL  CONTINUITY 

It  has  been  the  usual  practice  in  the  finite  element  solution  of  vector  field  problems  to  create  a  scalar 
approximation  for  each  vector  component  separately.  In  view  of  the  above  analysis,  we  must  do  some¬ 
thing  entirely  different:  we  must  generate  vector-valued  finite  elements  that  provide  continuity  of  only 
the  tangential  component  of  the  field.  No  finite  elements  having  this  property  have  been  published 
previously. 

In  the  two-dimensions,  the  idea  is  the  following:  consider  the  arbitrary  triangular  finite  element 

presented  in  Figure  1.  The  three  unit  vectors  I*.  ,  «  =  1,2,3,  parallel  to  the  three  sides  of  the  element 
are 

ri  -  Z  (*("„  +  «,T,)  (28) 

where  b.  and  c.  are  the  usual  parameters 

bi  -  Vn  -  y,2  (29) 


(31) 
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We  note  here  that  the  h.  may  also  be  written  in  terms  of  the  triangle  vertex  angles  9.  as 
h?  =  Aicotff^  +  cotOi2) 


where  A  is  plus  or  minus  twice  the  triangle  area 


A  =  c.„6.,  —  c.,6._ 

i2  il  il  i2 


(32) 


(33) 


The  dot  product  of  1  .  and  the  field  vector  "E  evaluated  on  side  *  is  the  tangential  component  of  ~E  on  that 
side.  It  is  the  quantity  that  we  wish  to  make  continuous  for  »  =  1,2,3. 

Now  note  that  the  tangential  component  of  ~E  evaluated  on  side  *  is  a  o"“-dimensional  scalar  and  can 
therefore  be  approximated  by  one-dimensional  N’th  order  interpolation  polynomials 

N 

r,-  •  ,  -  E  ‘X  (34) 

k=Q 

where  e'k  are  coefficients  and  0'k  are  N’th  order  interpolation  polynomials.  We  note  that  since  triangles 
have  three  sides,  and  since  there  are  (N+l)  parameters  on  each  side,  that  there  are  3(N+1)  coefficients 

i 

ek  ■ 

The  f)'k  may  be  written  in  terms  of  the  Silvester  polynomials 

P„M  -  i 


Pj!) 


m  ±  1 


(35) 


PN-k^il)  Pk^d 


(36) 


For  example,  if  t  =  1  and  N  =  2,  we  obtain 


Po  =  -  !) 

=  ?3(2?3  -  0  •  (37) 


The  location  of  the  interpolation  nodes  corresponding  to  these  polynomials  are  indicated  in  Figure  1(b). 
VECTOR  BASIS  FUNCTION 


Finite  elements  possessing  the  desired  tangential  continuity  are  provided  by  the  following  theorem: 


Theorem: 


Let  E  be  approximated  by  the  vector 


where 


N-l 


E 


(39) 


Then 


'tide  i 


N  ' 

-  E 


4 


(40) 


Proof: 


We  will  show  that  (40)  is  true  on  side  2;  the  proof  for  sides  2  and  3  is  very  similar.  Evaluating 
1  •  5  gives 


* 


E  =  2  {ClV*N-l  ~  C162V  ~  6lC3^-l  + 


1  ^2 


1  2  (  2  2  1 

2  \  {  "  C^N  + 


1  r  ,  3  3  ) 


(41) 


2  3 

On  side  1,  7^  =  0  and  =  0  by  construction.  Therefore  on  side  1  only  the  top  row  of  Equation 

(41)  survives.  However,  since  c^  -  6^  =  A  and  -  c^  +  6^  =  A,  Equation  (41)  yields 


1  1  *  E\tidt  1  ~<N- 1  +  7 


N 


(42) 


which  is  the  same  as  (40). 


THE  CASE  N  =  1 


For  N  =  1,  we  find  that  0Xg  =  ^  and  0*.  =  ^  .  Thus  (38)  becomes 


3  h. 


=  lz  L  2  {*^1  “ 


3  h. 


Fjr  E  2  {w,i  "  Wa} 


(43) 


Since  jf  —  —  (^  4-  bg  x  4-  C(.y),  we  can  rewrite  this  as 


,  * 

E  *  T  (/*  +  9*x  +  M  +  ~a  {fy  +V  +  \ 


(44) 
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Thus 


and 


fx  -  E  K6i2°il  -  C*i6,iai2) 

i=l 


g*  =  E  \  K  ~  ei6a6,2) 

1=1 


hT  =  E  A.  feo6,2C,l  -  ei6,lc,-2) 
»=1 


4  =  E  A.  (elc.2a,i  -  eic,xa,-2) 

i=l 


9%  =  E  \  -  elc,i6,2) 

i=i 


\  "  E  K  ~  elcac,'2) 

i=l 


V  •  E  =  -  +  hv) 


V  X  £  =  —  (g  -  h  )  1 

A  '  V  z'  z 


/  |V  X  E|!dl?  -  ij  (J>»  -  2„^  +  (,/) 


THE  CASE  N  =  0 


(45) 


(46) 


(47) 


(48) 


(49) 


Although  N  0  is  not  permitted  in  Equation  (38),  zeroth  order  tangential  continuous  finite  elements 
may  be  obtained  by  setting  e'g  =  e|  =  e'  in  Equation  (43).  The  result  is 
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1  3 

E  =  tE  *%  ~  *.i  «a) 


i«l 


1  3 

+  tE  e'\  (C,-2?,l  ~  C,l^) 


(50) 


i=i 


In  this  case  gx  =  0,  =  0  and  hx  =  —  gy=  Q  and  (44)  becomes 


1  1 

B  =  -J  (/,  +  »)  +  -J  Uy  ~  **) 


For  this  element  we  find  that 


V  •  E  =  0 


(51) 


V  X  "E  =  -  2g  z/A 


(52) 


(53) 


and 


/  IV  X  E|2dr?  =  i  ff2  .  (54) 

3-D  VECTOR  ELEMENT 

We  define  the  tangential  component  of  the  field  vector  at  side  t  to  be  positive  when  it  is  parallel  to  the 
unit  vector  /.  .  However,  to  have  a  consistent  sign  of  the  field  vector,  one  has  to  assign  the  tangential 
vectors  globally,  not  locally.  An  example  of  such  vector  mesh  is  shown  in  Figure  2.  It  is,  of  course, 
possible  to  have  other  vector  mesh  patterns  besides  the  one  shown;  these  patterns  are  related  through  a 
group  transformation.  For  example,  the  matrix  P  transforming  the  vector  (e1,  e2,  e3)  into  (-e1,  e2,  e3)  is 
defined  as 

/— 100  \ 

(e1,  e2,  e3)  P  =  (-e1,  e2,  e3)  where  P  =  [  010  ) 

\  001 J  (55) 


The  corresponding  element  matrices  on  the  transformed  mesh  are  given  as  follows: 
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Z .,  —  Z  —Z, 

jk  j  k 


' tjk  -  (4+4+4)5 


and  V  is  6  times  the  volume  of  the  tetrahedra.  Then 

'«  -  E  U,  «  -  £  (58> 

*=o 

where  1 . .  is  the  unit  edge  vector  on  the  edge  ij,  and 

«fxEl  1^=0  =  n.xE2  |£,.=0  (60) 

where  (1)  and  (2)  indicate  two  adjacent  tetrahedra  with  a  common  face  »",  and  n.  is  the  corresponding  unit 
normal  vector. 


THE  CASE  N  =  0 


In  this  case, 


Ex  =  yT.  7,.  lbtS~b,t)  *ij 
Fv  =  I'Z  • 


Let  us  define 


*  6  (C12’  *13’  *14’  *23’  *24’  *34^ 


Then,  the  matrix  elements  of  5  are 


(%,kl  =  3F  5,?n  ('i)  S‘>*  (*0 

+  yr  +  %  z*r)  • 


where  the  function  Sign  (ij)  =  1  when  i—jz i  2  ,  and  is  —1  otherwise,  and  ij  means  the  complementary 
index  of  ij e.g.,  12=34,  13=24,  and  34=12. 
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And  the  matrix  elements  of  T  are 


(tv . . ,  =  —  r-  7,.  (w7  /c.  -  w..  a:..  -  w.,  w,  +  w..  /c.) 

'  120V  'J  W  v  jl  tk  jk  1 1  1 l  jk  ik  p 


Where  W .  =  2  when  t  =  /,  1  otherwise,  and 

i)  J  ’ 


K..  =  bb.  +  c£.  +  d.d.  . 

>j  •  ]  >  J  •  J 


APPLICATION  TO  DIELECTRIC-LOADED  WAVEGUIDES 


L THEORY 


Assuming  the  wave  travels  in  Z  direction,  one  can  write  ET  as: 

E  =  +  j  ~ jAz 

where  e  t  =  +  e?  y  and  /9  is  the  propagation  constant. 

Substituting  (66)  into  (16)  yields 

F-  firr[?  'r.i2  +  'vx‘.i2  +  |VeJ2  -  ^  t*  +  %  S’ 


-<„*!  (F,l2  +  F/)  dn 


Defining 


'v  rp 

et  Tt  u 


tie 
z  z  ~z 


't  S  e 

t  t 


r  S  e 

z  z  ~z 

ej2  dn 


/  k  tl2  dr. 

/fj*  * 

/  I^x7(|2 
/|V 


(e~  ■  e~ )  H  (t:)  ”  hi •*  t,  +  ’*-£)in 


gives 
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F  =  FKc  -  kn  e' Me 

Ay  O  Av 


where 


=  (S'*  -  M 


K  =  }_  fp2Tt+St  (mA 
Hr  8W  st 


and 


M  = 


T*  % 

O  Tz 


The  first  derivative  with  respect  to  e“  gives 

K  e  —  if  M  e  . 


(70) 


(71) 


(72) 


(73) 


H.  RESULT 

Taking  a  rectangular  waveguide  half-filled  with  dielectric,  as  studied  by  Konrad,  as  an  example,  we 
approximate  the  transverse  component  of  the  field  vector  using  zeroth-order  vector  elements  and  the  lon¬ 
gitudinal  component  of  the  field  vector  using  first-order  scalar  element,  as  shown  in  Figure  1.  Unlike 
Konrad’s  result,  spurious  modes  do  not  occur  either  above  the  air  line  or  below  the  air.  It  is  interesting 
to  note  that  Konrad  had  incorrectly  identified  some  spurious  modes  as  physical  ones. 


APPLIED  TO  A  CAVITY  PROBLEM 


A  unit,  empty  cubic  cavity  was  modeled  by  the  zeroth-order  3-D  vector  element.  As  shown  in  Figure  ? 
the  agreement  is  reasonable.  The  result  also  shows  the  convergence  as  the  mesh  becomes  finer.  Opposed 
to  the  result  obtained  from  the  traditional  method,  the  convergence  is  simply  destroyed  by  the  fact  that 
the  number  of  spurious  modes  increases  with  the  number  of  the  degree  of  freedoms.  We  also  model  a 
dielectric-loaded  cavity  and  the  result  is  again  reasonably  good. 

THE  EVALUATION  OF  S*  AND  T* 


The  S  *  matrix  is  related  to  the  calculation  of 

I ivxei2  da  -  j  f 


where 


(74) 


9 


(75) 
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Let  us  define 

h  =  (hj  h2  h3)  and  e,  = 
Then, 


9  = 


e 


and 

2  ~  ~ 

g  =  e  h  . 

Since 


/i.  =  V(cot  on  +  cot  ei2)  , 


(76) 

(77) 


(78) 


(79) 


h  can  be  written  as 


h 


(80) 


and 

m  =  \/  cot  9  +  cot  6 ..  . 

•  il  «2 

Substituting  Equation  (81)  into  Equation  (80)  gives 


g2  —  A  e'  m  m  e 
From  Equation  (74),  S'  is 
S'  =2  mm. 


On  the  other  hand,  the  T'  matrix  is  related  to 


The  first  term  can  be  written  in  the  following  format: 


dn. 


(81) 


(82) 


(83) 


(84) 


(85) 
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3 

E  oA'iAi  - 

L,=i 


where 


K.  -I 


mimi^3^2*^1^3^  mim2^3^2*^1^3^  mim2^3^2*^2^3^ 

x(^1^3'^3^ 

tn2m2(b1?3‘b3$i)  m2m3(bl^3'b3^l) 

*(blW2) 


The  integral  of  is 

mimi(2b3-2b3Vf2b2> 

Ik,  «>  -  s  1  "'#">1(-2b,V2$ 

! 

[m3mi(*2b3br2b2) 

By  the  same  token,  the  second  term  becomes 


m2m2(2bj-2b1b3+b3) 

m3m2^"2b3b2’2bl)  m2m3(2bj-2b1b2+2b2) 


3 

E  VWa-'nta)?  “ 

‘■1=1 


and  the  integral  of  K2  is 


m1m1(2c3.2c3c2+2c2) 


f  K2  <tn  =  ^  m2m1(-2c1c2-2c3) 

ro3mi(-2c3cr2c2) 
By  applying  the  identities 


m2m2(2cf2ciC3+C3) 

m3m2(-C3C2-2c!) 


m3m3(2cr2ciC2+2c2) 


6-6.,+e.e.,  =  A 

i  il  i  il 

cot  $i2  i  5^  t1 


,2  2  . 

6,+C,  -  A 

(cot  9.^+cot  9i2) 


to  the  sum  of  and  ,  Eq.  becomes 
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24  C 


m1m1[2m1m1+6  cot  0j] 
m2m1[2  cot  03-2m3m3] 
m3m1(2  cot  6^- 2m2m2] 


therefore, 


A 

T* - 

24 


m2m2[2m2m2+6  cot  02] 
m3m2[2  cot  0^—  2m1m1] 


m3m3[2m3m3+6  cot  $3]  i 


(93) 

(94) 
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APPENDIX  C  -  C1  QUADRATIC  INTERPOLATION  FOR  DERIVATIVE 
CONTINUOUS  FINITE  ELEMENT  APPROXIMATIONS 

ABSTRACT 

A  highly  efficient  procedure  for  generating  smooth  surfaces  over  arbitrarily  spaced  data  points  is 
developed  in  this  appendix.  The  procedure  uses  quadratic  polynomials  for  the  construction  of  derivative 
continuous  surfaces  rather  than  the  cubic  polynomials  generally  employed  previously.  It  is  based  on  a 
subdivision  procedure,  dividing  each  triangle  in  a  triangulation  of  the  data  points  into  six  subtriangles 
and  fitting  a  quadratic  Bezier  surface  patch  over  each  subtriangle.  It  requires  only  function  and  first 
derivative  values  at  the  data  points  and  is  easily  implemented  by  means  of  simple  formulas  for  the  Bezier 
coefficients.  Since  two-dimensional  quadratic  polynomials  contain  only  six  terms,  while  ten  terms  are 
required  to  evaluate  a  cubic,  the  new  procedure  provides  a  major  improvement  in  the  efficiency  of  al¬ 
gorithms  for  representing  electromagnetic  fields. 
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INTRODUCTION 

The  generation  of  smooth  surfaces  that  interpolate  to  prescribed  values  of  a  two-dimensional  function  at 
an  arbitrary  point  set  is  one  of  the  central  topics  of  computer-aided  geometric  design  (CAGD).  Early 
work  focused  on  the  Coon’s  patch  as  the  preferred  method  of  generating  C J  surfaces  over  rectangular  and 
triangular  grids  I1);  more  recently,  because  of  their  easy  geometrical  interpretation,  techniques  for 
generating  smooth  surfaces  based  on  Bezier  polynomials  have  become  popular^. 

To  date,  however,  the  preponderance  of  work  in  this  area  has  been  based  on  the  considerable  flexibility 
of  cubics  for  providing  smooth  curves;  relatively  little  work  exists  for  using  quadratic  polynomials  to 
achieve  similar  results.  While  the  cubics  developed  previously  using  either  the  Coon’s  patch  or  the  Bezier 
technique  provide  excellent  results,  they  have  three  drawbacks: 

(1)  they  are  relatively  complicated; 

(2)  they  require  either  that  a  "twist  condition"  be  satisfied  at  the  data  points1  or  that  mid-side 
derivatives  be  given2;  and 

(3)  compared  to  quadratics,  they  are  expensive  to  compute. 

We  examine  here  the  requirements  for  C1  interpolation  by  quadratics  using  the  Bezier  technique  and 
show  that  a  simple,  efficient  and  accurate  algorithm  exists  for  generating  such  surfaces.  The  algorithm 
requires  only  function  and  derivative  values  at  an  arbitrary  set  of  points;  no  mid-side  or  twist  derviatives 
are  employed.  Further,  the  new  algorithm  is  considerably  faster  for  drawing  smooth  surfaces  than  are 
algorithms  based  on  evaluating  cubic  polynomials  since  fewer  operations  are  required  to  evaluate 
quadratic  polynomials  than  is  required  for  cubics.  Indeed,  it  may  be  argued  that  the  new  procedure  is 
optimal  for  C1  surface  interpolation  since  quadratics  are  the  lowest  order  polynomials  for  which  smooth 
interpolation  is  possible. 

The  pioneering  work  on  C1  interpolation  by  quadratic  polynomials  is  by  Powell  and  Sabin They 
showed  in  1977  that  C1  surfaces  could  be  produced  over  triangular  data  points  by  subdividing  the  original 
triangles  into  either  six  or  twelve  subtriangles.  Their  procedure  has,  however,  not  been  widely  employed 
due  to  its  complexity.  Sibson  and  Thomson4  have  also  introduced  a  quadratic  C1  interpolation  procedure 
but  their  method  is  limited  to  treating  only  rectangular  arrays  of  data  points.  The  method  developed  in 
this  paper  does  not  have  any  geometrical  restrictions  and  therefore  is  ideal  for  design  purposes.  The 
potential  of  quadratic  surface  patches  in  CAGD  has  also  been  recognized  and  is  discussed  in  a  recent 
paper  by  Sederberg  and  Anderson6. 

BEZIER  POLYNOMIALS 

The  Berastein-Besier  approach  to  surface  modeling  is  well  known  and  is  described  in  some  detail  in  a 
comprehensive  review  paper  by  Bohm,  Farin  and  Kahmann^l  With  regard  to  surface  modeling,  Bezier 
polynomials  provide  the  important  property  that  C1  continuity  can  be  ensured  between  triangular  patches 
by  requiring  that  each  pair  of  adjacent  subtriangles  in  the  Bezier  net  be  coplanar.  This  is  both 
conceptually  and  computationally  a  simple  task,  and  has  therefore  great  appeal  for  constructing  methods 
of  smooth  interpolation. 

For  quadratics,  the  Bernstein-Bezier  surface  patch  over  triangles  is  defined  as  follows 

/(*,v)  *  aw3)  -  E  s?  w 

ij,k  >  0 
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where  the  f.  are  the  homogeneous  or  barycentric  or  triangle  area  coordinates  of  projective  geometry^ 
and  the  coefficients  are  the  control  vertices  of  the  Bezier  net.  Note  that  unlike  the  case  with  ordinary 
Lagrange  interpolation,  the  coefficients  do  not  necessarily  provide  values  of  f(x,  y)  at  the  nodal 

points;  rather,  it  can  be  shown  that  the  surface  patch  will  lie  in  the  convex  hull  of  its  Bezier  net^8l 

The  constraints  on  the  control  vertices  in  neighboring  triangles  to  ensure  derivative  continuity 
across  the  interface  are  derived  in  Reference  2.  The  geometric  interpretation  of  these  constraints  is 
illustrated  in  Figure  1.  In  order  to  have  C1  continuity  between  quadratic  Bezier  patches,  each  pair  of 
subtriangles  of  the  Bezier  net  that  shares  a  common  boundary  —  the  two  shaded  quadrilaterals  in  the 
figure  —  must  be  planar. 

SUBDIVISION  PROCEDURE 

Consider  an  arbitrary  set  of  points  P  on  a  plane,  as  illustrated  in  the  example  in  Figure  2(a).  This 
point  set  may  be  triangulated  in  a  number  of  ways;  we  have  employed  Delaunay  triangulatioJ8'  to 
illustrate  a  possible  subdivision  in  Figure  2(b),  however  other  triangulations  may  be  employed  if  desired. 

We  assume  that  position  and  tangent  plane  data  are  given  on  P.  Let  Q  represent  a  new  set  of  points 
chosen  arbitrarily  so  that  each  triangle  T.  contains  one  and  only  one  point  Q.,  i*=l,  •  ■  •  ,  n(  where  n( 
represents  the  number  of  triangles  in  the  mesh.  By  connecting  each  point  Q.  with  the  associated  triangle 
vertices  and  with  the  three  neighboring  triangle  interior  points  Q.,  each  interior  triangle  is  subdivided 
into  six  triangles  as  shown  in  Figure  2(c).  Boundary  triangles  are  subdivided  into  six  by  adding  one  point 
B.  at  an  arbitrary  location  along  each  exterior  triangle  boundary. 

Consider  a  typical  pair  of  interior  triangles  as  illustrated  in  Figure  3.  The  quadratic  Bezier  net  further 
subdivides  these  triangles  as  indicated.  According  to  the  requirements  for  derivative  continuity,  the 
Bezier  net  of  each  of  the  shaded  areas  must  be  planar  in  order  to  achieve  C1  continuity.  The  question  is: 
Is  it  possible  to  construct  such  a  figure,  keeping  in  mind  the  interconnections  between  these  and  other 
elements? 

With  regard  to  the  original  triangle  vertices,  it  will  be  noticed  that  the  plane  of  the  Bezier  net 
surrounding  each  vertex  is  completely  defined  by  the  data  specified  at  the  vertex.  This  leaves  three 
shaded  regions  to  examine:  Two  triangular  patches  centered  on  the  interior  points  and  Qr  and  a 
diamond-shaped  region  efgh  spanning  the  interface  between  the  elements. 

The  triangular  patches  centered  on  and  Qt  are  easily  made  planar  by  requiring  th its  three  edges 
form  straight  lines  and  that  the  points  Q.  lie  in  the  plane  of  these  lines.  The  conditions  for  the  region 
efgh  at  the  center  of  the  common  triangle  edges  to  be  planar  are  more  complicated  and  are  provided  by 
the  following  lemma: 

Lemma.  Let  A,  B,  C  and  D  be  any  four  points  in  three-dimensional  space  as  shown  in  Figure  4. 
Further,  let  e,  /,  g  and  h  be  points  along  the  line  segments  AB,  Be,  CD  and  T)A,  respectively,  having 
the  properties  that  the  points  e  and  g  divide  the  line  segments  AS  and  DC  in  equal  ratios. 
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At  Dg 

a5  W 


(2) 


and  that  /  and  h  divide  BD  and  AD  in  equal  ratios. 


BT 

21 

5? 

aD 

(3) 


Then  the  four  points  t,  f,  g  and  h  are  coplanar. 

Proof:  Without  loss  of  generality,  we  may  adopt  a  Cartesian  coordinate  system  such  that  A  is  at  the 
origin,  B  lies  on  the  x  axis  and  C  lies  in  the  (x,y)  plane.  Then  the  coordinates  of  A,  B,  C  and  D  are 


A  =  (0,0,0) 

B  =  (b,0,0) 

C  «  (cvce0) 

D  -  (dvdrij  (4) 


Vectors  to  the  points  e-h  are  given  by 
7  —  TfjB  —  ilb~x 

7  —  +  (1-72)B 

-  (vi+i-V)”*+V2~i  • 

7  - 

-  (Vl+Wl)”.  +  (•71e8+d2-71rf2)ry 

+(d3-lld3)Tt 

r  «  7j25  -  ViW^r,  (5) 

We  need  to  show  that  the  vector  e  g  is  a  linear  combination  of  the  vectors  ejand  e7T 
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*  9  “  oe"7  +  0  7T  (6) 

The  x,  y  and  x  components  of  equation  (6)  are 

Vi+drVrV  —  °(vi+6-V-'yi6) 

+^(Vi~7i6) 

Va+Wg  -  o  Yj+^Vj 

Ws  =^2d3  (7) 

Rearranging  yields 

(7j-a72)Ci  +  (l-Tfj-^Jdj  -f  (-^-o+o^+a^+^Jfc  «=  o 

(71-072)c2  +  =  0 

(l“71~^72)d3  -  0  (8) 

Since  the  cf.  and  rf(  are  arbitrary,  it  follows  that 
7X  -  a72 

1  -  7j+072  (9) 

Solving  this  for  a  and  /9  gives 
°  -  7,/72 

^  “  (1— 7g)/72  (10) 

Since  b  is  also  arbitrary,  we  must  also  have  that 

-'/1-a+Of'T2+o'ri+^7l  —  0  (11) 


The  reader  may  verify  that  equation  (10)  also  satisfies  equation  (11). 
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By  construction,  the  quadrilateral  ABCD  satisfies  the  requirements  of  the  above  lemma.  It  follows  that 
C1  continuity  is  achieved  for  any  location  of  the  points  Q.  by  requiring  the  shaded  regions  around  the 
original  points  P.  and  around  the  interior  points  Qi  to  be  planar,  and  by  requiring  the  Bezier  node  E  to 
lie  in  the  plane  of  the  region  e/gh. 

A  FORMULA  FOR  C1  QUADRATICS 

We  now  proceed  to  derive  a  formula  for  the  Bezier  coefficients  4>i  on  the  subdivided  grid.  For  this 
purpose,  consider  the  triangle  T  shown  in  Figure  5.  In  this  figure,  Q  is  the  interior  point  in  triangle  T 
and  U,  V  and  W  are  the  interior  points  in  the  neighboring  triangles.  The  points  Q,  U,  V  and  W  may  be 
expressed  in  terms  of  the  Barycentric  coordinates  of  triangle  T  as 

Q  ■  (<Iv92,<]3) 

U  :  (t^.tij.ttj) 

V  :  (vvvvv3) 

W  :  (u;1(u>2,u»3).  (12) 

Since  the  points  U,  V  and  W  are  outside  triangle  T,  their  coordinates  will  not  be  limited  to  the  usual 
homogeneous  range  0<c<  1- 


LOCATING  THE  POINT  COORDINATES 

Our  first  task  is  to  determine  the  coordinates  of  the  Bezier  nodes  B..  To  simplify  the  algebra,  let  a,  0 
and  7  be  the  following  ratios 


3,7 


6,4 

*6,3 


where  d.j  represents  the  distance  between  nodes  i  and  j.  Now  note  that  the  line  Jbj  lj  is  described  by 
?,  —  («“?,)<  +  “  I’2'3 

for  0<<<1 .  At  node  B8,  fj  — =  0  so  that 

Thus,  substituting  (15)  into  (14)  gives  the  coordinates  of  as 
B6:  (0,  a,  1— o) 


(13) 

V 

(U) 

(15) 

(16) 


where 
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Vf«l“2 

o  ■=  . 

«1-«1 

In  the  same  way,  we  find  that 

B2:  (1-/J,  0,  0) 

B4-  H  l — y,  0) 

where 

0 


07) 


(18) 


V2“V3 


V"*2 


wi~h 


(19) 


The  coordinates  of  the  remaining  nodes  are  found  by  noting  that  they  are  located  midway  between 
nodes  with  known  coordinates.  For  example,  the  coordinates  of  £g  are  the  average  of  the  coordinates  of 


Bi  “d  Ba 
B* 


/1+T  1 — )r  \ 


2  ’  2  ’ 

*15  is  located  at  the  average  of  the  coordinates  of  B1  and 

B1&-  ((l+?1)/2,  9j/2,  q3/ 2j 

and  B,.  is  at  the  average  of  B.  and  B, 

' 

(  2  ’  2  ’  ?J 

Similar  results  are  obtained  for  the  other  nodal  coordinates. 
COMPUTING  FUNCTION  VALUES 


(20) 

(21) 

(22) 


By  definition,  we  are  given  the  function  values  and  the  values  of  its  x  and  y  derivatives  at  the  three 
vertices.  The  value  at  node  9  is 


^8  “  ^3  +  (^Ijt  1  3, 8^3,8 

where  1  .j  is  the  unit  vector  from  node  i  to  node  j 

Tn  -  t.  [(v*J  "»  + 

From  Equation  (20)  we  know  that 

*•  '  (t)*s  +  (rh 

».  -  (t2)*  +  (t2)*.  • 


(23) 


(24) 


(25) 
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It  follows  that 


(l — r) 


3-«  U. 


~  [(^s— X3^  “x  +  (»6-»s)  ~SI  ■ 
,9 


Thus,  Equation  (23)  yields 

*9  -  *^-rtxrxs)  f“l3  +  ^-7)(ys-v3)  ^la¬ 

in  the  same  way  we  find  that 

a  3d  8  di> 

*8  -  *3  +  2<I7'I3)  fet3  +  2  (V7"V3)  3^3 


(I— a)  v  3d, 

*6  +  "5“  t*7”*6)  ^5  + 


2  v'7  '3'  dy'  3 
5^l  ,  (1-a)  ,  \ 

A,'s  +  «  ^7  V&)  nJ* 


Tf  3d  •»  3d 

*10  “  *5  +  2  ^3~X6)  “d$h  +  2  ^3~y&^  ^5 

,  .  ^  1-0  ,  M,  A  (1-0  ,  *  «d. 

*13  “  *7  +  2  ^“^W7  +  2  1,3  1/7  3y 7 

a  3d  a  3d 

♦l2  -  *J  +  j  <'S-Z7fe  +  j  l*S-*7> 


1  3d  l  3d 

*15  “  *3  +  5  [(*l_1)x3  +  *2*5  +  «3X7l  ^3  +  2  [(,r1)y3  +  ^5  +  93y7]  ~l3 

1  3d  1  3d 

*17  -  *6  +  2  t*»*3  +  <«rlW  +  *3*yl  ^|6  +  -  l?,y3  +  (92-i)y6  +  93y7)  ^l5 

1  3d  1  3d 

*19  *  *7  +  2  (,»*3  +  *2*5  +  (93'1)*7J  ^7  +  &&  +  *2*5  +  ~|7 


The  remaining  nodal  values  are  computed  by  using  linear  interpolation.  First,  note  that 
*12,8  *19,18  *7,8 

*12,11  *19,17  *7,6 

Therefore,  the  line  segments  BnJB^,Bxi  and  B17rBl8.2?19  will  be  straight  provided  that 
*8  “  +  °*11 

*18  "  (l-Or)dw  4-  od17.  (31) 

Similarly 

*2  ■"  (l-0)d8  +  ^13 

*14  “  (l"Wl6  +  ^*19 

“  (1 — x)^J0  +  7d9 

*16  “  (i — r)*17  +  idj5-  (32) 

The  only  remaining  unknown  value  is  that  of  the  interior  node  Q,  alternatively  labeled  as  By  The 
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Bezier  net  of  this  node  must  lie  in  the  plane  of  the  triangle  E1&,£J7,£ig.  We  note  that  the  homogeneous 
coordinates  of  Q  are  identical  in  both  the  original  triangle  £3>£&,£7  and  in  the  central  subtriangle 
£1s,£17,£j9-  Therefore  the  value  of  Q  is 

*1  “  Mis  +  Ml7  +  Mio’  (33) 

In  view  of  the  above  lemma  and  the  proceeding  discussion,  we  have  established  the  following  theorem: 

Theorem:  Let  J  be  an  arbitrary  triangulation  of  an  arbitrary  set  of  points  £  in  a  plane.  Further,  let 
function  and  first  x  and  y  derivative  values  be  specified  on  P.  Then  a  piecewise  quadratic  Cl  surface  is 
formed  by  placing  a  point  Q.  at  an  arbitrary  location  inside  each  triangle,  thereby  subdividing  each 
triangle  into  six  triangles,  and  evaluating  the  Bezier  coefficients  of  the  quadratic  polynomials  on  the 
subdivided  triangles  according  to  the  formulas  in  equations  (27)-(29)  and  (31)-(33). 

SPECIAL  CASES 

While  the  above  theorem  provides  complete  flexibility  in  choosing  the  locations  of  the  interior 
subdivision  points  Q.,  it  is  appropriate  to  choose  specific  locations  for  these  points  in  computer 
implementations.  Two  logical  locations  for  these  points  are  (1)  the  triangle  centroids,  and  (2)  the  centers 
of  the  triangle  inscribed  circles.  We  note  that  although  the  centers  of  the  triangle  circumcircles  have  a 
special  meaning  in  the  commonly  used  Delaunay  triangulation  procedure,  these  circumcenters  represent  a 
poor  choice  for  our  algorithm  since  they  may  fall  outside  the  triangle  being  subdivided. 

THE  CENTROID  ALGORITHM 


If  Qi  is  chosen  to  be  the  centroid  of  triangle  t  then  its  homogeneous  coordinates  are 
1 
3 


92  *=  93  **=  -  In  this  case,  we  find  that 


uru2 

3^-1 


P 


V^3 

Sty-i 


W3-Wl 

3u»3-1 


(34) 


A  typical  surface  generated  by  the  centroid  algorithm  is  presented  in  Figure  6.  Figure  6(a)  shows  an 
irregular  triangulation  of  18  data  points  with  each  original  triangle  subdivided  into  six  triangles  by  using 
the  centroids.  Figure  6(b)  presents  a  typical  C1  quadratic  surface  produced  over  the  irregular  triangular 
grid  obtained  by  specifying  function  and  derivative  values  at  the  vertices.  As  is  apparent  from  this 
figure,  very  complicated  shapes  may  be  produced  by  using  only  a  few  elements. 


THE  INSCRIBED  CIRCUMCENTER  ALGORITHM 
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A s  is  well  known,  the  lines  connecting  the  center  of  an  inscribed  circle  to  the  triangle  vertices  bisects  the 
enclosed  angles  at  the  vertices.  Referring  to  Figure  7  and  using  the  area  properties  of  triangle 
homogeneous  coordinates,  we  And  that 


Aml  +  An2 


|A| 


(35) 


where  A.  is  the  area  of  the  subtriangle  indicated  in  Figure  7,  A  is  the  determinate 


V 3 


and  the  subscripts  m,  ml,  mS  are  cyclic  modulo  three.  Since  A. 

-  fV*  *ml  + 

Using  the  cotangent  identity 


2eot  6.,  Equation  (35)  yields 


eot  eml  +  cot  emi 


results  in  the  equation 


/2|A| 


1 


/|A|S 


(36) 

i 

(37) 

(38) 

(39) 


Equation  (39)  may  then  be  used  to  evaluate  a,  0,  and  in  Equations  (17)  and  (19). 


A  QUADRATIC  HONEYCOMB  ALGORITHM 

Far  in  has  published  a  C1  surface  generation  procedure  for  data  points  arranged  in  a  hexagonal  or 
"honeycomb"  pattern. ^  He  developed  this  procedure  by  using  the  properties  of  cubic  Besier  polynomials. 
We  shall  now  develop  a  similar  algorithm  based  on  the  properties  of  quadratic  Besier  polynomials. 

With  a  regular  array  of  hexagons,  a  triangulation  can  be  associated  as  described  in  Figure  8.  The 
vertices  of  the  hexagons  may  be  partitioned  into  the  two  groups  Cy  Cy  C&  and  Cy  Cy  Cy  We  can  then 
create  two  new  triangula«.ions,  both  of  which  consist  of  triangles  of  identical  site  and  shape,  one 
triangulation  based  on  the  points  Cy  Cy  C&  and  the  other  on  the  points  Cy  Cy  C#.  Following  the 
approach  of  Far  in,  we  let  each  group  of  hexagon  vertices  that  define  the  new  triangulations  form  a  plane 
in  three  dimensions  such  that  the  vertex  values  correspond  to  the  speciAed  values  on  the  honeycomb  grid 
By  insisting  that  the  control  vertices  of  the  besier  net  at  the  corners  of  the  original  triangle  lie  on  the 
average  of  these  two  planes,  we  can  express  each  control  vertex  as  a  weighted  sum  of  the  values  at  the 
honeycomb  vertices.  In  particular,  choosing  the  interior  points  Q.  to  be  the  triangle  centroids,  the 
following  simple  formulas  are  obtained 

*r  “  ;  <<?,+c1+cJ+c1+ci+cy 

*u  -  Jj  (iCj+sfCj+C,)-; 

*a  -  z  (nc,+c,)+4(c>+cs)+cs+c<  ). 


(«0) 
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By  symmetry,  identical  formulas  are  obtained  for  the  other  two  corners  of  the  original  triangle.  The 
remaining  control  vertices  are  obtained  by  enforcing  the  interior  C1  conditions  described  in  previous  sec¬ 
tions.  Thus  we  see  that  Farin’s  construction  is  equivalent  to  giving  position  and  tangent  plane  data  at 
the  triangle  vertices,  although,  in  fact,  no  derivative  values  are  specified  on  the  honeycomb  at  all. 

Figure  9  provides  an  example  of  surface  generated  by  the  quadratic  honeycomb  algorithm  where  only 
function  values  at  the  vertices  of  the  honeycomb  grid  were  specified. 

CONCLUSIONS 

Although  quadratics  have  often  been  dismissed  as  candidates  for  C1  data  interpolation,  they  are  in 
fact  ideal  for  such  application.  The  numerical  values  of  a  quadratic  C1  surface  are  controlled  purely  by 
function  and  first  derivative  values  at  the  data  points;  this  surface  is  defined  independently  of  the  mid¬ 
side  and  twist  derivative  conditions  required  by  cubics.  In  the  event  that  the  data  points  are  arranged  in 
a  regular  hexagonal  array,  a  very  simple  formula  exists  for  the  surface  purely  in  terms  of  surface  values  at 
the  honeycomb  vertices. 

The  main  advantage  of  the  procedure  developed  in  this  appendix  is  its  efficiency.  Quadratic 
polynomials  in  two-dimensions  contain  only  six  terms  compared  to  the  ten  terms  required  for  two- 
dimensional  cubics.  Although  the  algorithm  in  this  paper  requires  that  each  data  triangle  be  divided  into 
six  quadratic  subtriangles  to  ensure  C1  continuity,  while  previous  approaches  are  based  on  a  subdivision 
into  three  cubics,  the  number  of  subtriangles  required  by  the  surface  is  less  important  than  the  order  of 
the  polynomial  used.  The  reason  for  this  is  that  a  particular  subtriangle  in  a  surface  description  can  be 
located  ve.y  quickly;  a  much  greater  time  is  required  to  draw  the  surface  once  the  subtriangle  is  deter¬ 
mined. 

A  second  advantage  of  the  quadratic  approach  in  this  appendix  is  its  simplicity.  Surfaces  are  obtained 
by  generr  ing  an  arbitrary  triangulation  of  the  data  points,  subdividing  each  triangle  into  six,  and  using 
the  formulas  derived  in  this  paper.  Since  there  are  only  six  terms  in  a  quadratic,  the  formulas  defining 
the  quadr  itic  surfaces  defined  in  this  appendix  involve  relatively  few  fundamental  coefficients. 
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Figure  It  Two  neighboring  triangles  showing  the  associated  quadratic  Bexier 

nets  and  a  typical  control  surface.  The  shaded  areas  on  the  control  surface 
must  be  planar  for  derivative  continuity. 
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Figure  2:  (a)  A  representative  set  of  points  showing  the  Voronio  polygons. 

(b)  The  Delaunay  triangulation  of  the  points  in  (a),  (e)  Subdivision  of  the 
triangles  in  (b)  into  six  subtriangles. 
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Figure  6:  Node  numbering  scheme  for  *  triangle  subdivided  into  six  Bezier 

quadratics.  In  subtriangle  i,  the  node  numbers  are  given  by  the  formula 
{1,  (13+*),  (14+*  modB),  (l+i),  (7+i),  (2+*  mod6)} 


8:  (a)  The  triangulation  and  corresponding  Bezier  net  used  to  generate 

an  arbitrary  C'*  surface,  (b)  A  complex  C1  quadratic 
surface  obtained  by  specifying  random  function  and  derivative  values  at 

data  points. 
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APPENDIX  D  -  ELIMINATING  SPURIOUS  SOLUTIONS  IN  ELECTRO¬ 
MAGNETIC  WAVE  PROBLEMS  BY  USING  C1  FINITE  ELEMENTS 

ABSTRACT 

A  novel  approach  to  eliminating  spurious  or  nonphysical  solutions  in  finite  element  electromagnetic 
field  analysis  is  presented  in  this  appendix.  In  this  method,  C1  or  first  derivative  continuous  finite  ele¬ 
ments  are  used  to  construct  vector  field  solutions  that  are  exactly  divergence-free  both  within  and  across 
the  edges  of  the  elements.  The  appearance  of  spurious  modes  in  high  frequency  analysis  is  generally  at¬ 
tributed  to  the  inability  of  the  finite  element  field  solutions  to  satisfy  the  divergence  and  boundary  con¬ 
ditions  required  by  Maxwell’s  equations.  By  using  finite  element  vector  field  basis  functions  that  satisfy 
the  divergence  and  derivative  continuity  conditions  exactly,  spurious  solutions  are  completely  eliminated. 
In  this  appendix,  this  method  is  successfully  applied  to  inhomogeneous  dielectric  waveguide  problems.  In 
addition,  a  new  C1  triangular  finite  element  based  on  quadratics,  the  simplest  of  all  possible  general  C1 
elements,  is  introduced  for  finite  element  analysis. 

I.  INTRODUCTION 

The  fir’te  element  analysis  of  waveguiding  problems  is  one  of  the  most  important  areas  in  computa¬ 
tional  high  frequency  electromagnetics.  Different  methods  for  analyzing  dielectric-loaded  waveguides  and 
cavities  have  been  proposed.  They  include  vectoral  finite  element  formulations  in  terms  of  the 
longitudinal-electric  E*  and  magnetic  H*  field  components  [l]-[2]  as  well  as  formulations  that  are  based  on 
all  three  components  of  the  field  vectors  [3]-(5].  From  the  beginning,  it  was  recognized  that  a  major 
difficulty  in  these  analyses  is  the  appearance  of  spurious  or  nonphysical  solutions  [l]-[5].  These  spurious 
modes  intermingle  with  the  real  physical  solutions  and  make  the  distinction  between  correct  and  incorrect 
solutions  very  difficult.  One  cause  of  spurious  modes  is  that  the  finite  element  vector  solutions  do  not 

satisfy  the  divergence-free  condition  V  •  7}  —  0  [4]-[6]  that  is  required  by  the  physics  of  the  problem. 
To  remedy  this  difficulty,  various  methods  that  approximately  impose  the  zero-divergence  constraint  in 
the  finite  element  analysis  have  been  proposed.  The  penalty  function  method  [6]-[9],  for  instance,  incor¬ 
porates  the  condition  V  •  77  —  0  to  be  satisfied  in  the  least  square  sense  in  the  variational  functional. 
However,  this  method  requires  a  suitable  choice  for  the  penalty  coefficient.  Recently  a  method  that  ap¬ 
proximately  enforces  the  zero-divergence  condition  via  a  Galerkin  formation  has  been  introduced  [10]. 

Still  other  methods  that  are  based  on  approximating  V  •  77  =  0  by  finite  difference  [11]  or  by  using 
vector  field  basis  functions  that  are  approximately  divergence-free  [12]  have  also  been  proposed. 

Focusing  only  on  procedures  that  set  the  divergence  of  7?  to  zero  inside  finite  elements  is,  however, 
insufficient.  The  exact  solutions  of  electromagnetic  field  problems  are  not  only  nondivergent  inside  each 
finite  element,  they  are  also  nondivergent  at  inter-element  boundaries.  This  means  that  a  method  of 
correctly  specifying  the  derivative  of  the  solution  at  the  element  edges  must  be  developed. 

In  this  paper,  a  very  general  approach  for  constructing  finite  element  vector  field  basis  functions  that 
satisfy  the  zero-divergence  condition  both  inside  finite  elements  and  across  element  edges  is  presented. 
Using  this  method,  the  problem  of  spurious  solutions  in  waveguide  problems  is  completely  eliminated. 
The  novelty  of  this  approach  is  in  the  application  of  C1  or  first  derivative  continuous  finite  elements  in 
electromagnetics.  Although  the  Clough-Tocher  C1  cubic  [13]  is  well  known  in  structural  analysis,  the  use 
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of  C1  elements  has  thus  far  been  limited  to  mechanical  engineering  applications.  However,  in  this  paper, 
we  shall  demonstrate  that  Cl  elements  can  be  used  to  construct  vector  fields  that  are  exactly  divergence- 
free,  providing  an  alternative  way  of  eliminating  the  long-standing  problem  of  spurious  modes  in  com¬ 
putational  high  frequency  electromagnetics.  Further,  we  introduce  a  new  set  of  C1  basis  functions  based 
on  quadratic  polynomials  [14]-  [15]  for  finite  element  work.  These  C1  quadratics  have  the  advantage  of 
being  lower  order  than  the  standard  Clough-Tocher  triangle. 

H.  FORMULATION 

Consider  a  dielectric  waveguide  with  arbitrary  cross-section  Q  in  the  x-y  plane.  With  the  time  har¬ 
monic  variation  ejwt  assumed  throughout,  the  two  Maxwell  curl  equations  are: 


VxT?  =  —  jwno  7? 

(i) 

VxT?  =  jwto[tr]  2? 

(2) 

and  [er]  denotes  the  relative  permittivity  tensor  of  the  dielectric.  We  assume  throughout  that  the  relative 
permeability  of  the  dielectric  is  unity.  From  (l)  and  (2)  we  obtain  the  vector  wave  equation 


Vx 


([c^vxT?) 


fc2  7?  =  o 


(3) 


where  A:  =w  ji0f0-  In  finite  element  analysis,  this  equation  is  replaced  by  the  extremization  of  the  func¬ 
tional  [3] 


F  = 


([^-‘vxB)  -  k2  Tt  .  7/ 


<m  . 


(4) 


At  the  stationary  point  of  the  functional  (4),  we  obtain  (3)  as  the  Euler-Lagrange  equation.  The  exact 

solution  7/  which  satisfies  (3)  should  in  theory  also  satisfy  V  •  7?  =  0.  This  can  be  seen  by  taking  the 
divergence  of  both  sides  of  (3).  However,  problems  arise  when  we  obtain  an  approximate  solution  that 

does  not  satisfy  (3)  exactly  and  hence  V  •  7?  0.  Solutions  that  do  not  even  approximately  satisfy 

V  •  7?  —  0  are  generally  identified  as  spurious  modes.  To  eliminate  these  spurious  solutions,  one  can 

limit  the  subspace  of  trial  functions  7?  in  (4)  to  those  that  are  exactly  solenoidal  in  addition  to  satisfying 
the  essential  boundary  conditions. 


To  obtain  trial  functions  7?  that  satisfy  V 
to  expand  7?  in  general  as 


7?  =  0  a  priori,  it  is  mathematically  straightforward 


7?  = 


(5) 
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where  /9  is  the  propagation  constant  in  the  longitudinal  direction  and  ipx,  and  rPl  are  scalar  functions 
defined  over  J?.  While  the  preceding  expansion  is  formally  valid,  it  requires  that  each  of  the  scalar  func¬ 
tion  ip  have  at  least  continuous  first  derivatives,  that  is,  to  be  C1  throughout  the  entire  domain  f2.  This 

is  because  even  though  V  •  7?  as  0  holds  within  the  element,  continuity  of  7?  “cross  element  inter¬ 
faces  is  required.  Due  to  the  differentiation  in  (5),  C1  continuity  of  ip  automatically  implies  continuity  of 

7?  throughout  the  entire  domain  17.  In  other  words,  the  integral  form  of  Gauss’  law 


ds 


0 


(6) 


is  identically  satisfied  for  an  arbitrary  surface  e  within  the  waveguide,  including  surfaces  that  traverse 
element  edges. 

One  can  also  appreciate  the  need  for  C 1  continuity  by  substituting  (5)  into  (4)  and  so  that  F  is  written 
in  terms  of  the  trial  function  ip 

F  =  jf  |(VxVx?f  •  ([fr]_1VxVX?)  -  Jfc2(Vx^)*  •  (VXtft  d!2.  (7) 

The  functional  (7)  contains  second  derivatives  of  the  scalar  function  ip;  in  order  for  such  a  functional  to 
be  integrable,  C1  continuity  for  ip  is  required  [13].  It  is  emphasized  that  standard  C'  Lagrangian 

polynomials  are  inadequate  to  yield  an  identically  zero-divergence  field  77;  such  fields  require  at  least  C1 
continuity. 

m.  QUADRATIC  C1  TRIANGULAR  FINITE  ELEMENTS 

The  analysis  presented  in  this  paper  is  based  on  the  recently  developed  C1  quadratic  triangular  finite 
element  that  is  described  in  [14]-[15] .  Its  construction  is  based  on  a  subdivision  principle.  An  interior 
point  Q,  such  as  the  centroid,  or  more  generally  the  triangle  incenter,  is  selected  for  each  element  in  the 
mesh.  Then  by  joining  each  interior  point  to  its  element  vertices  and  the  interior  points  of  neighboring 
elements,  each  element  is  subdivided  into  6  subtriangles  (Fig.  1).  Over  each  subtriangle,  a  second  order 
polynomial  is  defined  using  the  Bernstein-Bezier  representation  [16].  After  imposing  C 1  continuity,  the 
resulting  C1  element  interpolates  to  the  function  and  to  the  first  x  and  y  derivatives  at  the  element  ver¬ 
tices. 

For  the  purpose  of  this  paper,  it  is  sufficient  to  summarize  the  linear  relationship  between  the  19 
Bezier  coefficients  <p,  and  the  function  and  first  derivatives  data  of  <j>(x,y )  at  the  triangle  vertices  as  given 
in  Table  1. 

IV.  FINITE  ELEMENT  FORMULATION 

In  conventional  finite  element  analysis,  discretization  of  the  functional  (4)  results  in  the  general  eigen¬ 
value  problem  [3] 


(S]{//}  -  k2[J]{H) 


(8) 
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where  the  components  of  the  column  vector  {//}  are  the  values  of  H ^  H '  and  H 'g  at  the  nodal  points  and 
the  matrices  [S]  and  (7)  are  given  by 

{H}+[S\{H}  =  [  (VxTT*  .  dfi  (9) 

Jn 

{H}+[T]{H}  =  f  (75)*  .  (75)  dn.  (10) 

Jn 

According  to  (5),  the  vector  {H)  is  related  to  a  column  vector  {P}  which  contains  the  parameters 
describing  the  scalar  C1  functions  if>.  ,  *  =  x>y,z •  This  is  represented  by  the  “curl*  matrix  [C\ 

{H}  =  [C](P>  .  (11) 


Substituting  (ll)  into  (9)  and  (10)  results  in  a  modified  eigenvalue  problem  that  guarantees  V  •  75=0 


[cf[s](c]{p}  =  k2\atT[T][C]{P) 


(12) 


For  the  Cl  quadratic  described  in  the  previous  section,  the  entries  of  {P}  are  the  function  and  first 
derivatives  of  V  at  the  vertices  of  the  element.  The  matrix  [C\  is  constructed  for  each  element  using  (5) 
and  the  information  provided  in  Table  1.  The  modified  [5]  and  [7]  matrices  in  (12)  are  then  computed 
for  each  element  and  assembled  to  form  the  global  matrices  for  the  whole  problem. 

An  equivalent  procedure  is  to  carry  out  the  expansion  in  (7).  The  modified  [S']  and  [T]  matrices  can 
then  be  computed  directly  without  resorting  to  the  congruence  transformation  in  (12). 

The  extremization  of  the  functional  (4)  requires  the  essential  boundary  condition  75  •  a  —  0  to  be 

satisfied  along  perfectly  conducting  surfaces.  These  boundary  conditions  on  75  then  correspond  to  con¬ 
straints  on  the  V*’s  in  (&)•  Examples  will  be  given  in  the  next  sections  on  the  application  of  these  con¬ 
straints. 

V.  TEST  EXAMPLES 

A.  EMPTY  RECTANGULAR  WAVEGUIDES 

In  the  first  example,  we  shall  illustrate  the  present  method  by  solving  for  the  eigenmodes  of  an  empty 
square  metallic  waveguide  of  width  W.  For  the  TM  modes  of  the  empty  waveguide,  we  let  V’X==V'  ~0  and 
write 


75  =  Vx*,(*,y)e-**  7  . 


(13) 


The  essential  boundary  condition  75  •  a  =0  then  implies  that  V\(>  •  a.  =  0  where  a.  denotes 

*  it  t 

the  tangent  vector  to  the  metallic  boundary  Ton  the  x-y  plane.  To  enforce  this  condition,  it  is  sufficient 
to  impose  =  0  on  P.  In  the  same  way,  to  solve  for  the  TE  modes  we  write 
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~E  =  Vx^(i,y)e  j0‘  az  (14) 

and  require  that  E  •  a  (  =  Oonf.  This  is  imposed  by  setting  •  a  =  Oonf. 

The  waveguide  was  discretized  by  using  81  C1  elements;  the  first  10  TM  modes  computed  are  given  in 
Table  2.  It  is  interesting  to  note  that  even  with  a  homogeneous  guide,  the  classical  analysis  based  on  (8) 
yields  several  zero  eigenvalues  in  addition  to  the  correct  nonzero  values.  It  is  these  zero  eigenvalues  that 
shift  to  nonzero  spurious  frequencies  if  the  guide  is  partially  filled  with  dielectric.  In  contrast,  the  present 
method  as  based  on  (12)  does  not  produce  any  "zero"  eigenvalues  with  an  empty  waveguide.  Every 
eigenvalue  obtained  corresponds  to  a  physical  solution.  The  next  example  shows  that  even  with  in¬ 
homogeneous  dielectric  waveguides,  no  spurious  modes  appear  in  the  analysis. 

B.  DIELECTRIC-LOADED  WAVEGUIDES 

Consider  a  square  metallic  waveguide  half-filled  with  dielectric.  The  dielectric  has  relative  permitivity 

cr  =  1.5  and  the  geometry  is  shown  in  Fig.  2.  Solving  for  the  LSM  modes  of  the  guide  we  let 

=  rb  =  0  and  write 
if  z 

H  =  VXtl>z{x,y)e~j0z  a‘z  .  (15) 

The  boundary  conditions  7/  •  a  n  —  0  on  f  now  correspond  to  ij>z  —  0  along  y  —  0  and  y  —  W. 

In  one  test,  100  C1  elements  were  used  yielding  a  total  of  319  unknowns.  For  0W  =  10,  the  first  ten 
eigenvalues  from  this  test  are  given  in  Table  3  and  compared  to  those  obtained  by  Koshiba  et  al  [8],  All 
solutions  obtained  by  using  the  C1  procedure  correspond  to  physically  correct  solu„.  ms  so  that  the 
problem  of  spurious  modes  is  completely  eliminated. 

VI.  CONCLUSION 

A  novel  approach  for  eliminating  spurious  modes  in  waveguiding  problems  is  presented.  This  method 
is  based  on  the  use  of  C1  finite  elements  for  constructing  vector  field  solutions  that  are  divergence-free. 
Unlike  previous  approaches  that  impose  the  zero-divergence  condition  approximately  by  modifying  the 
variational  procedure,  field  solutions  from  the  present  method  satisfy  the  zero-divergence  constraint  ex¬ 
actly  both  inside  and  across  the  edges  of  the  finite  elements.  It  is  found  that  in  this  way,  spurious  solu¬ 
tions  do  not  occur. 

This  paper  also  introduces  the  use  of  C1  quadratic  triangles  to  finite  element  analysis.  These  quad¬ 
ratics  are  the  simplest  C1  element  possible  and  hence  the  most  economical  to  use  in  many  applications. 

It  should  be  noted  that  the  method  presented  in  this  paper  for  constructing  divergence-free  finite  ele¬ 
ment  vector  fields  is  a  very  general  one.  Its  application  is  not  limited  to  the  two-dimensional  waveguide 
problems  considered.  For  instance,  we  are  presently  investigating  the  elimination  of  spurious  modes  in 
3-D  electromagnetic  cavities  by  applying  three-dimensional  C1  finite  elements.  The  procedure  can  also  be 
used  to  generate  zero  curl  vector  fields  by  changing  the  curl  operator  in  (5)  to  the  gradient  operator. 
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Table  Is  Formulas  for  the  Bezier  Coefficients  of  Composite  C 1 
Quadratics  Finite  Elements 


4* 

9$ 

9  $ 

3<J 

9$ 

9$ 

«  j 

3y  1 
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^7 
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1 

«2 

b2 

$n 

a 

a*2 

ab2 

(l-<x) 

(1-a'iaj 

(l-a)b, 

$19 

1 

“3 

b, 

where  (q1(q2,q3)  denotes  the  simplex  coordinates  of  Q  and 


a2=l  K  V^V1)  x5+%  x7] 

-  a4=i  (X7"X3) 

•  a8=i  (X3_X5) 

’  a8=I  (X8_X7) 

°  =  Vd7,5 

^  ~  d3,2^d3,7 
7  =  d5,4/d5,3 

The  expressions  for  b(  are  obtained  from  those  of  a(.  by  substituting  y^.  for  x.  . 
d  ..  denotes  distance  between  node  t  and  node  j. 


al==2  x3+q2  xs+%  *7! 

a3=5  [qi  X3+V5+(%-^  X71 

U-7)  ,  . 

a5=—  (XS“X3) 

(»-•)  ,  , 

*7=—  (X7“Xs) 


(»-«  ,  V 

a9=—  (X3~X7) 
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Table  2:  The  First  Ten  TM  Modes  in  &  Square  Waveguide 


TM  Mode  C1  Finite  Element  Exact 


1,1 

4.4552 

4.4429 

1,2 

7.0873 

7.02481 

2,1 

7.0894 

7.02481 

2,2 

8.9847 

8.8858 

1,3 

10.1472 

9.93459 

3,1 

10.1475 

9.93459 

2,3 

11.5391 

11.3272 

3,2 

11  5751 

11.3272 

1,4 

13.4612 

12.9531 

4,1 

13.4672 

12.9531 

Table  3:  The  First  Ten  LSM  Modes  for  0W  =  10  in  the 
Dielectric  Waveguide  Shown  in  Figure  2. 


No. 

C1  Finite  Element 

Hayata  et  al. 

1 

8.8109 

8.8093 

2 

9.8834 

9.3896 

3 

10.2751 

10.2752 

4 

11.1073 

11.1038 

5 

11.3835 

11.2677 

6 

11.4685 

11.4501 

7 

12.2654 

11.9882 

8 

12.7957 

12.6686 

9 

12.9510 

12.8092 

10 

13.4298 

12.9575 
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APPENDIX  E  -  MODELING  THE  OPEN  BOUNDARY  CONDITION  IN  THE 
FINITE  ELEMENT  SOLUTION  OF  ELECTROMAGENTIC  SCATTERING 
PROBLEMS 

ABSTRACT 

Electromagnetic  scattering  problems  are  solved  in  this  appendix  by  using  the  finite  element  method 
within  an  artificial  bound  picture  frame,  and  analytical  expressions  in  the  exterior  region.  Two  schemes 
for  coupling  the  solutions  in  the  two  regions  are  proposed,  namely,  a  "two-boundary  approach"  and  a 
"transfinite  elements  approach."  The  two  boundary  approach,  as  implied  by  the  name,  employs  two 
boundaries  in  the  solution  procedure  and  the  superposition  theorem  to  determine  the  proper  linear  com¬ 
bination  of  modes  for  the  problem.  By  allowing  the  exterior  region  and  the  finite  element  solution  region 
to  overlap,  the  coefficients  in  the  modal  expansion  are  determined.  As  for  the  transfinite  element  method, 
it  is  based  on  variational  principles,  choosing  different  basis  functions  for  the  finite  elements  solution 
region  and  for  the  exterior  region.  The  application  of  the  two  boundary  approach  to  waveguide  junction 
problems  as  well  as  to  object  scattering  problems  is  presented  here.  Application  of  the  transfinite  element 
method  to  electrostatic  problems  and  to  object  scattering  problems  are  also  studies. 
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SECTION  1  INTRODUCTION 

1.1  OPEN  BOUNDARY  SCATTERING  PROBLEMS 

A  number  of  important  electromagnetic  problems  involve  the  bistatic  scattering  of 
electromagnetic  Helds  from  metal  and  dielectric  objects  in  open  regions.  The  primary  categories  of 
these  scattering  problems  are  the  following: 

1.  Radiation  from  antennas. 

2.  Electromagnetic  scattering  from  arbitrarily  shaped  and  inhomogeneous  penetrable 
bodies. 

3.  Microwave  planar  integrated  circuits. 

4.  Waveguide  junctions. 

In  these  problems,  knowledge  of  the  Held  distribution  is  important  in  two  different  regions:  (1)  at 
a  very  great  distance  from  the  scattering  object  and  (2)  within  the  vicinity  of  the  scatterer  itself. 

Since  the  geometry  of  electromagnetic  scatterers  is  in  general  complicated  and  often  three- 
dimensional,  numerical  methods  are  preferred  for  their  analysis.  These  methods  may  be  divided 
into  two  groups:  (1)  integral  equation  methods,  and  (2)  differentia]  equation  methods.  Of  these 
two  approaches,  by  far  the  most  commonly  used  has  been  the  integral  approach,  popularized  by 
Harrington  with  his  monograph  “Fields  Computation  by  Moment  Methods-.1  However,  despite 
their  popularity,  integral  methods  have  failed  to  solve  a  number  of  important  scattering 
problems,  most  notably  those  involving  large  or  complicated  variations  in  permeability  and/or 
permittivity. 

On  the  other  hand,  in  other  engineering  disciplines,  differential  methods  have  been  very 
successful  in  modeling  large,  complicated  inhomogeneous  structures.  Differential  methods  appear 
to  have  an  advantage  precisely  where  integral  methods  are  weak  -  namely  in  modeling  the 
interface  condition  between  media  with  different  properties.  They  have  had  limited  application  to 
scattering  problems  for  only  one  reason:  differential  methods  cannot  accomodate  easily  the  open 
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or  infinite  region  encountered  in  scattering  problems.  While  a  number  of  procedures  to  model 
open  boundary  conditions  with  differential  methods  have  been  proposed  and  will  be  reviewed 
below,  none  of  the  existing  methods  is  completely  satisfactory  in  terms  of  simplicity,  ease  of  use, 
and  efficiency. 

The  purpose  of  this  appendix  is  to  develop  new  methods  of  modeling  electromagnetic  scattering 
from  complicated  shapes.  We  employ  the  finite  element  method  to  model  the  fields  in  the  near 
field  scattering  region  and  develop  two  new  ways  to  model  the  open  region  enclosing  the  scatterer. 
This  appendix  is  restricted  to  considering  electromagnetic  scattering  problems  to  two  dimensions. 
We  specifically  treat  the  scattering  that  appears  at  waveguide  junctions  and  from  cylindrical 
scatterers  in  free  space.  However,  the  methods  developed  in  this  report  will  also  be  useful  for 
solving  the  general  electromagnetic  scattering  problems. 

1.2.  MODELING  THE  OPEN  BOUNDARY  CONDITION 

Since  finite  element  methods  are  best  suited  to  modeling  electromagnetic  fields  in  the  immediate 
vicinity  of  the  scatterer,  some  method  of  representing  the  fields  in  the  infinite  region  must  be 
developed,  as  well  as  procedures  for  coupling  the  finite  element  solution  region  to  the  infinite 
region.  With  regard  to  waveguide  junction  scattering,  such  procedures  were  developed  relatively 
early  and  are  known  by  the  name  modal  analysis  technique23  The  idea  behind  this  method  is  to 
replace  the  waveguide  discontinuity  problem  by  two  appropriate  semi-infinite  uniform  waveguides 
separated  by  a  so-called  "discontinuity  region".  Muilwyk  and  Davies5  adapted  an  experimental 

a 

procedure  proposed  by  Montgomery  et  al  for  finding  the  scattering  parameters  of  a  discontinuity 
by  locating  the  nulls  of  its  standing  wave  patterns.  They  solved  the  field  problems  in  the 
discontinuity  region  by  using  the  finite-difference  method.  A  deficiency  in  the  method  proposed 
by  Muilwyk  and  Davies  is  that  finding  the  locations  of  the  nulls  in  the  standing  wave  pattern  is  a 
trial  and  error  procedure  and  on  each  trial  an  eigenvalue  problem  must  be  solved. 

Marin7  has  applied  the  finite  element  method  to  solve  the  scattering  problems  with  lossy 
scatterers.  A  nonlocal  boundary  condition  was  employed  to  derive  a  variational  formulation  of 
the  scattering  problem,  and  the  finite  element  method  was  applied  to  determine  the 
approximation  to  the  near  fields.  Scattering  amplitudes  are  determined  by  means  of  an  integral 
representation  obtained  from  Green's  formula  and  properties  of  the  nonlocal  boundary  operator. 
But  the  complexity  of  the  calculation  for  the  nonlocal  boundary  operator  has  limited  the  usage  of 
this  method.  Recently  Koshiba  and  Suzuki8  proposed  a  better  way  of  using  the  integral  technique 
for  solving  the  problem  of  a  H-plane  waveguide  junction  containing  lossy  ferrite  poets  of  arbitrary 
shape.  The  procedure  is  to  consider  two  planes  in  each  port  and  then  to  relate  the  unknowns  on 
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these  two  planes  using  an  integral  equation  procedure.  The  result  is  a  set  of  equations  that  are 
combined  with  the  equations  obtained  from  the  finite  element  method  and  solved  directly  for  the 
field  values  inside  the  discontinuity  region.  Unfortunately  the  resulting  matrix  equation  is 
complex  and  non-symmetric  even  for  the  lossless  cases. 

With  regard  to  open  cylindrical  scattering  problems,  a  common  feature  »-  numerical  modeling 
methods  is  the  subdivision  of  the  problem  into  a  finite,  closed  region  called  a  "picture  frame*  and 
the  surrounding  open  space.  A  finite  element  or  finite  difference  approximation  is  then  made  for 
the  fields  within  the  picture  frame  and  some  device  used  to  couple  the  exterior  field  boundary 
condition  to  the  picture  frame  boundary. 

In  one  approach,  proposed  by  Silvester  and  Hsieh®  for  electrostatics  problems,  a  single  picture 
frame  is  defined  and  the  energy  functional  of  the  entire  region  exterior  to  the  picture  frame 
evaluated  and  added  to  the  interior  region  energy  functional.  The  solution  is  therefore  involves 
the  electrostatic  field  in  all  space,  but  by  a  judicious  choice  of  exterior  field  approximation 
functions  explicitly  solves  for  the  field  only  in  the  region  interior  to  the  picture  frame.  A  major 
difficulty  with  this  formulation,  however  is  that  the  integrand  along  the  picture  frame  boundary 
is  singular  and  difficult  to  compute.  In  another  approach,  proposed  by  McDonald  and  Wexler10, 
two  concentric  picture  frames  are  defined  and  the  integral  equation  relating  the  potentials 
between  the  two  picture  frames  is  used  to  specify  the  boundary  conditions.  Fields  outside  of  the 
outer  picture  frame  are  never  considered  in  the  solution  process;  the  integral  equation  n  erely 
provides  a  relationship  between  internal  field  values.  However,  the  energy  functional  in  the 
exterior  region  can  be  evaluated  by  using  an  integral  transformation  and  weighted  Gaussian 
quadrature  formulas.  The  programming  requirements  of  this  procedure  are  also  relatively  difficult 
and  this  has  limited  its  application.  By  using  the  mode  decomposition  technique,  Chang  and  Mei 
11  presented  a  method  which  they  called  "unimoment  method"  to  calculate  the  scattered  fields  of 
dielectric  cylinders  of  arbitrary  cross  section  or  of  inhomogeneous  material.  The  basic  technique  of 
the  method  is  to  use  the  finite  element  method  inside  a  circular  picture  frame,  and  to  expand  the 
fields  outside  in  cylindrical  harmonics.  Coupling  th'  ‘r'vior  and  exterior  solutions  is  performed 
by  matching  the  interior  and  exterior  fields  with  rc‘  •'<  to  both  the  function  values  and  the 
normal  derivative  values  on  the  circle.  There  are  two  disadvantages  with  this  method:  first  for  N 
expansion  functions  in  the  exterior  region,  2*  N  boundary  value  problems  must  be  solved  as  the 
basis  functions;  second,  derivatives  of  the  finite  element  solutions  for  each  mode  on  the  boundary 
must  be  calculated. 

Bettess12  in  1977  proposed  a  procedure  he  dubbed  "infinite  elements".  In  this  method,  a  series 
of  shape  functions  analogous  to  Lagrange  polynomials  but  including  an  exponential  decay  term 
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are  generated.  The  approximate  interpolation  functions  are  unbounded  geometrically  hence  the 
name  "infinite  elements*.  Unfortunately  the  accuracy  of  this  method  depends  on  the  ability  of 
the  user  to  place  the  right  decay  factor  in  the  infinite  elements.  Also  forming  such  infinite 
elements  is  relatively  inefficient  and  complex  so  that  usage  of  this  method  is  relatively  limited. 

A  somewhat  better  group  of  methods,  are  called  "ballooning"  and  its  variants.  Ballooning  was 
invented  by  Silvester  for  modeling  the  Laplacian  problems  of  the  open  type.  It  also  employs  a 
picture  frame  to  limit  the  size  of  the  finite  element  mesh;  however,  in  this  case  the  exterior 
boundary  is  modeled  by  using  uniform  layers  of  finite  elements  on  the  outside.  Unfortunately, 
this  method  has  deficiencies  as  well:  first  the  algorithm  is  not  efficient  when  applied  to  the  wave 
equation;  second  the  field  values  outside  the  picture  frame  can  not  be  obtained  easily.  Dasgupta 
14  extended  the  ballooning  algorithm  into  the  cloning  algorithm.  In  cloning,  the  formulation  has 
led  to  a  quadratic  matrix  equation.  The  unknown  relates  to  the  desired  matrix  pertaining  to  the 
infinite  or  semi-infinite  region.  The  major  contribution  of  the  cloning  method  is  that  no 
analytical  expressions  are  required  to  supply  the  Sommerfeld  radiation  condition.  But  in  order  to 
obtain  the  unknown  matrix,  a  set  of  quadratic  eigenvalue  problems  has  to  be  solved.  More 
recently,  a  method  called  "infinitesimal  scaling"  has  been  proposed  by  Hurwitz15  16  17.  The 
infinitesimal  scaling  method  seeks  to  determine  the  coefficient  matrix  which  expresses  the 
contribution  of  an  element  to  the  global  functional  in  terms  of  the  solution  values  on  the 
boundary  nodes  by  using  a  scaling  concept  analogous  to  that  used  in  ballooning,  but  in  the  limit 
that  the  size  of  the  boundary  layer  element  approaches  zero.  By  taking  this  limit,  a  nonlinear 
first-order  differential  equation  for  the  coefficient  matrix  is  obtained.  The  price  that  is  paid  is 
the  need  to  solve  the  a  nonlinear  matrix  differential  equation. 

1.3.  THE  NEW  PICTURE  FRAME  PROCEDURES 

In  this  appendix  we  develop  two  new  picture  frame  techniques.  However,  we  limit  the  shape 
of  the  boundary  between  the  picture  frame  and  the  open  region  to  be  simple;  as  circle  in  the 
cylindrical  scattering  problems  and  lines  in  the  waveguide  junctions.  By  using  the  trir  jgular 
finite  element  method  to  solve  the  fields  interior  to  the  simple  shaped  picture  frame,  a  simple, 
efficient  and  stable  algorithm  for  numerical  calculations  is  produced.  We  shall  call  the  two  new 
coupling  schemes  described  in  this  report  the  "Two-boundary"  approach  and  the  "Transfinite 
elements"  approach. 


The  idea  behind  the  "Two- boundary"  approach  is  derived  from  the  modal  analysis  technique  in 
which  the  finite  element  method  and  the  superposition  theorem  are  used  to  calculate  the  fields. 
First  the  modal  analysis  technique  is  used  to  get  the  analytical  expressions  for  the  fields  in  the 
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region  exterior  to  the  picture  frame.  Second,  a  number  of  boundary  value  problems  in  the  picture 
frame  are  set  up  and  solved  by  the  finite  element  method.  In  this  step,  finite  element  solutions 
provide  an  expression  for  the  numerical  solutions.  Finally,  we  require  that  the  analytical  solution 
and  the  numerical  solution  match  along  the  second  boundary.  The  detailed  description  of  the 
■Two- boundary"  approach  will  be  covered  in  Section  3.  It’s  application  to  the  waveguide 
junction  problems  as  well  as  the  results  are  in  Section  4,  and  its  application  to  open  2-D 
scattering  problems  is  given  in  Section  5. 

For  the  •Transfinite  elements,"  we  construct  basis  functions  outside  the  picture  frame  in  such  a 
way  so  as  to  convert  the  integral  of  the  functional  for  the  whole  region  into  a  integral  over  just 
the  picture  frame  plus  a  boundary  integral.  Applying  the  variational  principle  then  provides  a 
matrix  equation  to  be  solved  for  the  scattering  solution.  Section  6  presents  the  application  of  the 
"Transfinite  elements"  method  to  model  simple  electrostatic  problems.  The  application  to  the  2-D 
scattering  problems  is  investigated  in  Section  7. 
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SECTION  2  INTRODUCTION  TO  THE  FINITE  ELEMENT  METHOD 


The  basic  concept  of  finite  element  analysis  consists  of  three  parts  : 

1.  Interpolation:  Distributed  physical  qualities  may  be  approximated  by  using  linear 
combination  of  interpolating  functions. 

2.  Mesh  generation:  A  complicated  physical  shape  may  be  subdivided  into  many 
subregions,  and  the  integrals  of  distributed  functions  over  the  whole  region  may  be 
broken  into  summations  of  the  integral  over  the  subregion. 

3.  Variational  principles:  For  each  problem,  derive  a  corresponding  functional  either  by 
using  Galerkin’s  method  or  by  using  the  Rayleigh-Ritz  method.  Then  apply  the 
variational  procedure  to  the  functional  to  derive  a  matrix  equation. 


2.1.  VARIATIONAL  PRINCIPLES 

The  finite  element  method  is  based  on  variational  principles  developed  for  the  most  part  in  the 
last  century.  Variational  principles  for  electromagnetic  field  problems  can  be  derived  from  two 
alternative  points  of  view,  namely,  the  Galerkin’s  method  and  the  Rayleigh-Ritz  method. 

Consider  the  following  two  dimensional  Helmholtz  equation 
V2  rl>  +  k2 =  p(x,y)  in  /? 

rl>  —  0  on  T  (2.1) 

involving  what  may  be  assumed  to  be  a  scalar  potential  variable  The  quantity  k2  is  a  constant 
invariant  with  position,  whilst  p  is  a  given  forcing  function  over  the  problem  region  Cl,  and  r  is 
the  boundary  of  Cl.  For  simplicity,  we  only  consider  the  case  of  Dirchlet  boundary  conditions 
here. 
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2.1.1.  GALERKIN’S  METHOD 

1  ft 

The  steps  in  Galerkin’s  method  to  solve  (2.1)  are  as  follows  : 

1.  Find  a  complete  function  space  Au  to  be  called  the  “trial  function  space*.  Any 
function  u  €  Au  must  satisfy  the  Dirchlet  boundary  conditions  on  F.  The  trial 
function  space  A  will  be  used  to  write  the  numerical  approximation  of  i/>  of  (2.1)  in 
the  form 

N 

»-=l 

where  N  is  the  dimension  of  A  ,  the  tl>.’ s  are  the  unknown  coefficients  and  a.(x,y)’s  are 

0  1  1 

the  basis  functions  of  A  . 

2.  For  every  trial  function  u  in  the  trial  function  space  Au,  define  a  residual  function 

/^(x.y)  —  V~u(x  ,y)  +  *2  u(x,y)  -  />(  x,y)  (2.3) 


3.  Construct  another  function  space  Ay  to  be  called  the  “testing  function  space".  Any 
function  in  Ay  will  be  called  a  testing  function  and  will  be  used  to  weight  the  residual 
function  fu(x,y)  defined  in  (2.3).  Although,  doesn’t  have  to  be  the  same  as  Au, 
often  it  is  convenient  to  have  Aa  m  Ay . 

4.  Assign  to  each  pair  (v,u),  v  €  Ay  and  u  €  a  scalar  number  according  to  the  bilinear 
form  B(v,u) 

B{v,u)  =  -Jn  t <z,y)  ft(x, y)  <K2  (2.4) 

By  Green’s  theorem  B(v,u)  can  be  written  as 

B(v,u)  *=  V  wV  u  (ffi  —  k2  Jn  vu  <fJ7  4-  v  p  dfi  —  JpV‘fa<*r  (2-5) 

Because  for  any  testing  function  v  €  Ay,  v  =  0  on  F,  this  simplies  to 

JB(v,u)  =  Jn  V  v*  V  u  dfi—  ^  Jp  vu  ^  +  Jn  v  P  dfi  (2.6) 

5.  The  problem  of  finding  the  solution  of  (2.1)  is  thus  converted  into  the  following 

Find  €  A%  such  that 

B{v,*)  =  0  *  for  V  v  6  (2.7) 

This  solution  is  called  a  “weak  solution"  because  the  solution  of  (2.7)  will  be  the 
solution  of  (2.1)  only  in  an  approximate  sense. 
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6.  Write 

N  M 

£  *iai*>v)  v==Yl  vj&lx,y) 

•“1  J»“l 

where  M  is  the  dimension  of  Ay,  and  the  #.(x,y)’s  are  basis  functions  of  vi^.  Then 

E  E  «><  jnv)i‘v‘<-k,0fi  *.+E  *7  /„  V1"7 


(2.8) 


j“  1  i««l 

7.  Because  =  0  for  V  v  g  dy,  we  have 

N 


r=l 


e  { /0K*v“<-*v<)‘fDh+ jjifdn-0 

j= 


(2.9) 


i=l 


8.  Usually,  du  m  /1v,  then  (2.9)  can  be  written  in  the  matrix  form 

( s-  jfc2r)  -  7  (2.10) 

where 

is  the  column  vector  of  the  coefficients  ^ 

i 

S,  Tare  matrices  with  5..=  [  Va.»Va.df} 

•j  Jn  j  • 

T. .  =  f  a  a .  dfl 

13  Jn  J  • 

and  j  is  &  column  vector  with  /.  =  —  /  a. pdf} 

Jn  'y 

Galerkin’s  method  is  analyzed  in  detail  in  reference18. 

2.1.2.  RAYLEIGH-RITZ  METHOD 

The  Rayleigh-Ritz  method  is  more  complicated  than  Galerkin’s  method  because  it  is  not  derived 
directly  from  the  differential  equation.  Rather,  one  first  forms  a  functional  from  the  stored 
energy  in  the  system;  minimizing  the  functional  over  the  set  of  approximating  functions  then 
gives  the  approximate  solution.  The  concepts  embedded  in  the  Rayleigh-Ritz  method  can  be 
explained  in  the  following  steps. 

1.  The  energy  functional  F(^)  corresponds  to  equation  (2.1)  is 
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<V,p>  +  <p,'P>  -  <!f',VV>  -k2<<P,4'> 

=  [  *pdn+  f  p*dn-  f  vv^dn- k2 [  wdn 

jo  Jo"  Jn  Jn 


(2.11) 


where  the  inner  product  <u,v>  is  defined  as 


<u,v>  =  I  u  v 

Jo 


<m 


Again,  by  Green’s  theorem  and  applying  the  same  boundary  conditions  as  in  the 
Galer kin’s  method 


Jn*pdn+  Jnp*<m+  jnvt*V'i'dn-k2  (2.12) 

2.  Since  we  are  interested  in  variations  in  F(tf')  about  the  solution  point,  we  may  begin 
by  writing  the  arbitrary  function  ^  as  a  sum  of  the  true  solution  rp  and  a  second 
arbitrary  function  £  multiplied  by  a  scalar  e,  ^  =  ip  +  (£ 

The  first  variation  of  F(tf')  about  the  solution  ip  is  given  by 

«WI^  =  flt=0 


which  gives 


=  2  Jntpdn+ 2  Jnvt.V4>dn-k2Jnwdn 

by  Green’s  theorem 


=  -2  Jnt(v24,  +  k2ii,-p)dn 

Since  ip  is  the  solution  of  equation  (2.1),  the  first  variation  of  F(if')  about  the  solution 
is  zero.  Thus  the  energy  functional  F(^)  is  stationary  about  the  true  solution  ip. 

3.  If  the  energy  functional  F(tf')  is  stationary  about  a  function  <p,  i.e. 

*  =  <f>  +  i  i 

«WU  =  fUo  =  ° 


which  gives 


=  -2  |oe(v2*  +  *2*-p)dr?  =  0 

Since  £  is  an  arbtrary  function,  we  have  V2  <p  +  k2<p  -  p  —  0.  So,  the  function  4>  will 
also  be  the  solution  of  (2.1). 
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4.  Approximating  by  a  linear  combination  of  interpolation  functions  or. 
N 

»=1 

gives  the  energy  functional  as 
N 


(2.13) 


w  =  2 5Z  jnaipdn+H  H  { jnVai*Var k2°iaj  ^ 
1=1  1=1  r=l 


5.  Setting  the  first  variation  of  F(lP)  to  zero  in  order  to  find  the  stationary  point,  then 
results  in 


E{  f  V  a  v  a  .  —  k2a .  a  .}  dSJ  'if .  +  f  a  p  dfl  —  0 

In  1  i  »  y  j  Jn  • 

i=l 

i  =  1....N 

This  matrix  equation  then  is  the  same  as  the  one  in  (2.10) 

The  details  of  the  Rayleigh-Ritz  method  can  be  found  in  reference19. 

2.2.  MESH  GENERATION 

The  finite  element  method  involves  breaking  the  problem  region  into  mesh  of  small  subregions. 
Each  subregion  is  called  an  element. 

The  mesh  can  made  in  any  number  of  ways.  One  common  method  is  to  use  a  square  or 
rectangular  grid  with  uniform  sized  elements.  The  method  that  has  been  used  in  this  thesis 
employs  triangular  elements  of  varying  sizes.  One  advantage  of  triangular  grids  is  their  great 
flexibility  in  fitting  the  geometry  of  different  problems.  Also  when  combined  with  the  Lagrange 
interpolation  polynomials,  to  be  discussed  in  the  next  section,  the  integration  over  the  elements 
can  be  performed  analytically10.  The  algorithm  used  to  generate  the  triangular  mesh  is  the 
Delaunay  triangulation  procedure  developed  earlier  at  CMU20.  Use  of  the  Delaunay  algorithm 
guarantees  that  the  triangulation  has  a  maximum  sum  of  the  smallest  angles  of  the  triangles. 
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2.2.1  DELAUNAY  ALGORITHM 

For  an  arbitrary  set  of  points,  a  geometric  structure  called  Voronoi  polygons  is  defined  by  the 
regions  where  every  point  within  a  polygon  is  closer  to  that  point  than  to  any  other.  The 
intersection  of  three  or  more  Voronoi  polygons  is  called  a  vertex.  This  vertex  is  by  definition 
equidistant  from  the  three  or  more  original  points  associated  with  the  polygons  that  formed  the 
vertex.  The  Delaunay  triangulation  of  the  system  is  formed  by  forming  triangles  from  the  three 
original  points  associated  with  each  vertex. 

One  property  of  the  Delaunay  triangulation  is  that  the  circumcircle  of  a  Delaunay  triangle  may 
not  contain  a  fourth  Delaunay  vertex.  This  property  is  the  basis  of  the  swapping  algorithm21  for 
Delaunay  triangulation  which  leads  itself  readily  to  finite  element  applications.  An  arbitrary 
triangulation  is  formed  from  the  original  points.  For  each  triangle,  the  neighboring  point 
associated  with  each  of  the  adjacent  triangles  is  checked.  If  any  of  those  points  is  within  the 
circumcircle  of  the  first  triangle,  the  diagonal  is  swapped.  See  Figure  2-1. 


Figure  2-1:  Swapping  a  diagonal 

In  the  event  that  a  diagonal  is  swapped,  it  is  necessary  to  check  the  pew  triangles  again  for 
additional  swaps.  This  procedure  is  repeated  until  no  more  swaps  are  necessary. 

2.3.  INTERPOLATION  FUNCTIONS 

The  family  of  polynomials  is  one  of  the  simplest  class  of  functions  to  employ  as  interpolation 
functions  and  were  the  first  to  be  applied  in  the  finite  element  method.  In  this  thesis,  we  only 
consider  the  most  common  interpolation  polynomials  used  in  the  finite  element  method,  namely 
the  Lagrange  interpolation  polynomials. 

If  or.(x,y)  is  the  Lagrange  interpolation  polynomial  associated  with  point  (Xj,y4),  then  we  have  by 
definition 
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(0  if  j  , 

°  ^=1  ifj-i  (2‘M) 

Equation  (2.14)  leads  to  the  following  physical  interpretation  of  the  coefficients  in  (2.13) 

N 

*— i 

so  that  the  coefficient  if',  is  also  the  function  value  at  node  (x.,y  ). 

j  v  j’V 

In  this  thesis,  we  use  second  order  triangular  elements.  For  each  triangle,  there  are  six  nodes 
and  six  quadratic  Lagrange  interpolation  polynomials.  Inside  each  triangle,  we  interpolate  the 
function  by  these  six  interpolation  polynomials  as 
8 

*{x,y)  =  £  *,(*.y) 

1-1 

To  derive  the  interpolation  polynomials,  it  is  convenient  to  introduce  "homogeneous  coordinates" 
on  the  triangle. 

2.3.1.  HOMOGENEOUS  COORDINATES 

Consider  an  arbitrary  triangle  with  vertices  numbered  (1,2,3)  and  vertex  coordinates  (Xj.yj), 
(x2,y2),  (*3,y3)  “  indicated  in  Figure  2-2  ,  and  let  P  be  a  point  located  within  this  triangle. 

The  coordinates  (x,y)  of  the  point  P  may  be  written  in  a  symmetric  form  with  respect  to  the 
triangle  (1,2,3)  by  defining  a  new  type  of  coordinate  system  called  homogeneous  coordinates  in 
the  triangle.  These  homogeneous  coordinates  are  defined  as  follows:  Let  dj  be  the  distance  of  the 
point  P  from  the  base  2-3  of  the  triangle,  h}  be  the  corresponding  altitude,  and  ^  the  ratio 


The  quantity  ^  is  called  a  homogeneous  coordinate  for  the  point  P  because  its  range  is  from  0  to 
1  for  any  point  within  the  triangle  (1,2,3). 

In  the  same  way,  homogeneous  coordinates  f2  and  f3  ranging  from  0  on  the  other  two  sides  to  1 
at  vertices  2  and  3  are  defined  by  the  equation 

4. 

C  =  — -  ,  m  =  1,2,3 

m  h 

m 

In  this  equation,  dm  and  hm  are  defined  analogously  to  d{  and  hj. 


(*»,  w 


Figure  2-2:  A  typical  triangle  with  vertices  numbered  (1,2,3).  dj  is  the 

distance  of  the  point  P  from  the  base  and  hj  is  the  triangle  altitude. 

The  triangle  (1,2,3)  is  divided  into  three  by  the  point  P. 

The  homogeneous  coordinates  f  are  also  called  triangle  area  cordinates  Tor  the  following 
reason:  The  area  of  a  triangle  is  one-half  its  base  times  its  height.  Accordingly,  the  ratio  of  the 
area  Ap23  of  the  triangle  (P,2,3)  to  the  area  Aj  „3  of  the  triangle  (1,2,3)  is 


Consequently,  the  homogeneous  coordinate  ^  represents  not  only  the  ratio  of  the  distance  d}  to 
the  altitude  hj,  but  also  the  ratio  of  the  area  of  the  shaded  triangle  (P,2,3)  in  Figure  2-2  to  the 
area  of  the  original  triangle  (1,2,3). 

The  three  homogeneous  coordinates  fj,  f3  cannot  be  independent,  since  only  two  linearly 
independent  coordinates  can  exist  on  a  plane.  The  relationship  between  the  (fm)  is  revealed  by 
examining  the  properties  of  the  triangles  in  Figure  2-2  where  (1,2,3)  represents  the  original 
triangle  and  the  remaining  three  triangles  are  composed  of  the  lines  connecting  the  point  P  with 
vertices  1,  2,  and  3.  The  area  A]U  of  the  original  triangle 

Apto  +  Alp3  +  Al2p  =  ^123 

Dividing  both  side  of  this  equation  by  AJ23  yields 

fl+f2+?3=1 

Thus,  the  sum  of  the  triangle  area  coordinates  is  equal  to  unity.33 
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2.3.2.  INTERPOLATION  ON  TRIANGLES 

Now,  consider  the  problem  of  defining  quadratic  interpolation  polynomials  on  the  triangles. 
Since  three  point  values  are  required  to  define  a  quadratic  polynomial,  quadratic  interpolation 
polynomials  must  have  three  interpolation  nodes  on  each  side,  see  Figure  2-3. 


Figure  2-3:  Node  ordering  for  a  second  order  triangular  element 


The  polynomial 


Qj(x,y)  is  xero  along  the  lines  =  0  and  =  1/2  and  it  is  one  at  fj.  It  follows  that  aj(x,y) 
equals  xero  at  the  five  points  (f^  f2>  f3)  =  (1/2,  1/2,  0),  (1/2,  0,  1/2),  (0,  1,  0),  (0,  1/2,  1/2),  (0, 
0,  1)  and  is  equal  to  one  at  the  vertex  point  number  1.  Second-order  polynomials  which 
interpolate  to  1  at  the  other  two  triangle  vertices  are  defined  similarly. 


To  form  a  polynomial  which  interpolates  to  one  at  the  center  of  side  1-2,  we  need  to  define  a 
polynomial  which  is  zero  along  side  1-3  and  2-3  and  has  a  unit  value  at  the  point  (1/2,  1/2,  0). 
The  polynomial 

a4(*,y)  = 

provides  these  properties.  Polynomials  which  interpolate  at  the  remaining  two  midside  nodes  are 
defined  similarly.22 
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SECTION  3  THE  TWO-BOUNDARY  APPROACH 


As  is  implied  by  the  name,  two-boundary  method  for  solving  open  Held  problems  requires  two 
surfaces  to  be  defined  in  the  solution  process.  Because  two  boundary  surfaces  are  used,  space  is 
divided  into  three  parts  by  this  algorithm.  Figure  3-1  shows  a  portion  of  the  infinite  x-y  plane  in 
Cartesian  coordinates.  A  bounded  region  called  R}  is  defined  so  as  to  contain  all  sources, 
inhomogeneities,  and  anisotropies.  The  boundary  of  is  called  Fj,  and  another  larger  boundary 
Fj  is  placed  around  Fj.  The  region  in  between  Fj  and  F2  is  called  R2,  and  the  region  exterior  to 
F2  is  called  R^.  The  region  enclosed  by  F 2,  i.e.  Rj  ©  R2,  will  be  the  region  where  the  finite 
element  method  is  used,  which  is  then  called  the  'finite  element  solution  region". 


Figure  3-1:  Separation  of  the  problem  region 

8.1.  MODAL  ANALYSIS 

Since  the  region  exterior  to  Rj,  i.e.  R2  ©  R?;  is  homogeneous  and  source  free,  homogeneous 
solutions  of  the  Helmholtz  equation  can  be  found  by  using  the  method  of  separation  of  variables 
technique.  In  the  region  R2  ©  R(,  the  analytical  expression  for  the  fields  can  be  approximated 
by  M  modes  as 

M 

i-l 


(3.1) 
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where  a.’s  are  the  coefficients  that  need  to  be  determined  and  the  7-(x,y)  represent  the  analytical 
solution  for  mode  i  in  the  system. 

Since  equation  (3.1)  is  valid  for  the  region  R2  ®  Re,  it  also  holds  for  the  fields  on  the  boundary' 

rr 

M 

*lr2  =  £  a.  VI,s')|r2  M 

t—i 

This  leads  to  the  following  boundary  value  problem  for  the  fields  in  Rj  ®  R2 

(  V2  +  Jfc2 )  #  —  0  in  ®  /?2 

M 

]r  ai^x>y)  i  r  on  r2  (3-3) 

i— 1  2 

Equation  (3.3)  is  not  a  solvable  boundary  value  problem,  because  the  boundary  conditions  depend 
on  the  coefficients  a.'s  which  are  not  known  yet. 

We  may  use  the  superposition  theorem  to  decompose  equation  (3.3)  into  M  boundary  value 
problems  corresponding  to  the  modes  7.(x,y).  The  fields  pattern  for  each  mode  7.  is  designated  A. 
and  can  be  obtained  by  using  the  standard  finite  element  method  to  solve  for  the  fields  in  region 
Rj  ®  R2  with  the  boundary  conditions  A.(x,y)  =  7.(x,y)|^  on  T’j. 

Once  all  the  A.(x,y)’s  are  determined,  the  fields  inside  the  region  Rj  ®  R2  can  be  expressed  in 
terms  of  the  unknown  coefficients  a.’s 

l 

M 

*(*.v)  =  x/z,y)  f°r  iz<v )  €  z?j  ©  z?2  (3.4) 

i«=l 


8.2.  MATCHING 

Although,  we  have  not  yet  solved  the  fields,  the  only  problem  remaining  is  to  determine  the  M 
coefficients  a/s. 

Notice  that  we  have  two  expressions  for  the  solution,  one  is  equation  (3.1)  which  is  valid  in 
region  Rj  ®  R(,  and  the  other  is  equation  (3.4)  which  is  valid  in  region  Rj  ®  Rj.  It  follows  that 
both  solutions  are  valid  in  the  common  region  R2  and  that  they  should  have  same  values  on  the 
boundary  7^. 

Matching  the  modes  obtained  from  these  two  solutions  on  the  boundary  7^,  either  by  the 
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weighted  residue  method 
solved  for  the  coefficients 
can  be  calculated  either 
calculations  are  presented 


or  the  least  square  method,  results  in  a  matrix  equation  which  can  be 
a's.  Once  the  coefficients  a/s  are  known,  the  fields  values  at  any  point 
by  using  equation  (3.1)  or  by  using  equation  (3.4).  Details  of  the 
in  Sections  4  and  5. 
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SECTION  4  TWO  BOUNDARY  APPROACH  IN  RECTANGULAR 
WAVEGUIDE  DISCONTINUITIES 


4.1.  THE  NATURE  OF  THE  PROBLEM 


Figure  4-1:  A  general  two-port  rectangular  waveguide  discontinuity 

A  microwave  network  with  an  arbitrary  discontinuity  and  two  rectangular  waveguide  ports  is 
pictured  in  Figure  4-1.  The  discontinuity  may  be  an  arbitrary  discontinuity  in  a  single  guide  or 
the  coupling  between  two  different  guides.  Many  filters,  matching  sections,  phase-correction  units, 
and  other  components  can  be  represented  by  this  structure. 

Figure  4-1  consists  of  three  parts:  rectangular  waveguide  I  on  the  left,  rectangular  waveguide  II 
on  the  right  and  the  region  between  them  called  the  discontinuity  region  (D.R.).  Note  that  guides 
I  and  II  need  not  be  colinear  and  the  discontinuity  region  can  be  of  arbitrary  shape  and  may 
contain  dielectric  objects  or  even  metal  objects;  however  the  overall  structure  must  be  uniform 
(ie.  have  constant  cross  section)  in  one  direction  (either  the  "broad*  or  "narrow"  transverse 
direction).  This  requirement  simplifies  the  resulting  boundary  value  problem  to  two-dimensions. 

If  a  transverse  plane  wave  is  incident  from  left,  then  there  will  be  a  reflected  wave  in  guide  I 
and  a  transmitted  wave  in  guide  II  and  evanescent  modes  in  both.  The  polarization  of  the  waves 
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in  the  structure  provides  two  fundment&l  kinds  of  waveguide  discontinuities:  H-plane 
discontinuities  in  which  the  H  field  is  in  the  plane  of  the  discontinuity  and  E-plane  discontinuities 
where  the  E  field  is  in  the  plane  of  the  discontinuity. 

4.2.  FORMULATION  OF  THE  TWO  BOUNDARY  APPROACH 


zl  z3  z4  z2 


Figure  4-2:  Finite  Element  Solution  Region 

Assume  that  the  discontinuity  in  Figure  4-1  is  on  the  x-*  plane  and  that  the  dimensions  of  the 
guides  are  such  that  only  the  dominant  mode  can  propagate  in  each  guide.  Because  the 
formulation  is  similar  for  H-plane  discontinuities  and  E-plane  discontinuities,  only  the  H-plane 
discontinuity  formulation  is  presented  here.  In  Figure  4-2,  we  have  placed  the  boundary  planes 
z=zl  and  z=z2  to  form  the  "finite  element  solution  region”  as  discussed  in  Chapter  3.  The 
finite  element  solution  region  must  enclose  the  discontinuity  region  in  the  waveguide.  Also  shown 
in  Figure  4-2  are  planes  z>=z3  and  z=z4  which  served  as  the  matching  boundaries.  For  an  H- 
plane  discontinuity,  with  an  incident  wave  HJ0  from  the  left  ,  the  incident  electric  field  Ejnc  has 
only  a  y  component  and  can  be  expressed  as 

Ein  c=«'n(H«p{-y^).  /=v42-H2  (4.1) 

where  k  is  the  wavenumber  and  ^  is  the  propagation  constant  in  guide  I. 

The  modes  for  the  region  exterior  to  the  finite  element  solution  region  are: 
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j  tnx  j 

t(*.*)  =  *«'"( - )  •  «*p(o,  *) 

a 


»n  guide  I. 


II  tnx  II 

TT ,.  (*,*)-=•*«(— )«ezp(a|.  z)  in  guide  II. 

where 

^  >I/2  “d 

i  =  1... 


(4.2) 


The  modes  tt.1  represent  the  reflected  wave  modes  in  the  waveguide  I,  and  Tr.n  are  the  transmitted 
wave  modes  in  waveguide  II.  Set 


tf(x)  =  ««"(—)  rf*)  =  exPi  a!*) 

and 

V-ft*)  =  «'n(~)  ^\z)  —  exp{  aj7z ) 


(4.3) 


Then  the  modes  can  be  written  as  the  product  of  two  functions,  one  depends  only  on  x  and  the 
other  depends  only  on  z. 


tJ(x,z)  =  V$x)  •  ^.{z) 
and 

Tj(x,z)  —  tf-7(x)  •  ^(z) 


(4.4) 


Combining  equations  (4.1)  —  (4.3)  gives  the  analytical  expression  for  the  fields  inside  guides  I  and 
II  as 

.  N 

E*(x,z)  =  V»{(*)  •  {  exp{-jl3rz)  +  R^x{z)  j  +  ^  R^j(x)^j(z)  z  <  z^ 


and 


i'—2 


&'{*,*)  -  £  '(*>»?(*> 


(4.5) 


i-l 


Here  the  R.  and  T  are  the  reflection  and  transmission  coefficients  for  each  mode,  N  is  the  number 
of  modes  used  in  guide  I  and  M  is  the  number  of  modes  used  in  guide  II. 

From  the  discussion  in  Chapter  3,  two  set  of  fields  patterns  \}  and  X.D  within  the  finite  element 
solution  region  are  obtained.  Both  solve  the  Helmholtz  equation  in  the  finite  element  solution 
region  ,  but  with  different  boundary  conditions  on  planes  z  *=  Zj,  z  =  Zj.  The  boundary 
conditions  for  X.1  are  X^x.Zj)  ^(x)  »nd  X.^x.Zj)  *  0,  and  for  X.n,  X^x.Zj)  «=  0  and 
X.°(x,z2)  *—  V.n(x).  It  follows  that  the  numerical  solution  in  the  finite  element  solution  region  may 
be  written  as 


113 


N  M 

-  r« jK-isfy+ff/u,) |x(+£  Rfaj  x'  +  £  rVV2)  x"  (4.6) 

'  '  «'— 2  i-l 

Now,  we  introduce  two  inner  products  for  complex  valued  functions  one  on  the  plsne  t  — -  i3, 
the  other  on  the  plsne  »  —  i<t  m 

<«(*).«(*)>/  “  J0  v’(x)u(x)dx 

snd  (4.7) 

rD  . 

<  “(*  W*)  >  //  “  J0  v(x)u(x)dx 

where  represents  the  complex  conjugate.  Since  both  equations  (4.5)  snd  (4.6)  hold  in  the  region 
*5  >  t  >  ij  snd  *2  >  t  >  »4,  which  then  give  the  following  equations  need  to  be  solved  for  the 
coefficients  R.  snd  T.. 

<  TZ  (*,x3) ,  d>'(x)  >!  —  <  ^(x.Xj) ,  V>f(x)  >t  i  **=  1....N 

snd 

<  'E(x,xi),i>'\x)>ll— <En(x,xi),i/>[\x)>n  »  —  (4.8) 


Finally,  substitute  equation  (4.5)  snd  (4.6)  into  (4.8)  results  in  the  (N+M)  times  (N+M)  complex 
matrix  equation 
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4.3.  NUMERICAL  RESULTS  AND  COMPARISON 

A  number  of  waveguide  discontinuity  problems  have  been  studied  and  the  results  from  this 
analysis  will  be  compared  with  that  obtained  from  other  theories  and  from  experimental  data. 
The  results  here  employ  three  modes  on  both  sides  for  matching.  Our  experience  indicates  that 
evanescent  modes  die  out  very  quickly  so  that  even  if  a  small  finite  element  solution  region  is 
used  "three  modes  matching"  still  provides  good  results. 

4.3.1.  H-PLANE  STEP  AND  E-PLANE  STEP 

Figure  4-3  shows  a  rectangular  waveguide  of  width  "a"  joined  symmetrically  with  a  guide  of 
width  "D".  The  discontinuity  in  this  problem  occurs  in  the  H-plane.  Figure  4-3  provides  a 
comparison  of  our  results  with  the  approximate  analytical  solution  of  Lewin23.  For  the 
discontinuity  in  the  E-plane,  the  results  from  present  method  as  well  as  those  from  Marcuvitz2* 
are  given  in  Figure  4-4. 

4.3.2.  THICK  IRISES 

Here  we  consider  two  types  of  irises  in  rectangular  waveguides.  Figure  4-5(a)  and  Figure  4-5(b) 
show  the  magnitude  and  phase  angles  of  the  reflection  coefficient  due  to  square  inductive  irises  of 
various  sizes.  The  agreement  between  the  results  computed  here  and  that  obtained  by  Davies5  is 
within  a  few  percent.  Figure  4-6(a)  and  Figure  4-6(b)  show  the  real  and  imaginary  parts  of  the  E 
field  contour  lines  for  6/ a  =  0.2,  a/X  =  0.7. 

Figure  4-7  compares  the  VSWR  ratio  calculated  from  present  method  with  the  analytic 
solutions  from  Kerns25  for  the  single  half-round  inductive  obstacle  for  r/a  =  0.2. 

4.3.3.  WAVEGUIDE  BENDINGS 

The  method  we  have  developed  may  also  be  used  to  calculate  the  scattering  of  waves  by 
waveguide  bends.  The  structure  of  the  waveguide  bends  we  are  considering  is  shown  in  Figure 
4-8.  In  this  figure,  the  planes  "t"  are  reference  planes  for  the  parameters  R  and  T  .  Table  4-1 
gives  a  comparison  of  present  theory  and  experimental  data  from  Marcuvitz2*  for  different 
bending  angles. 
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4.3.4.  DIELECTRIC  LOADED  WAVEGUIDE 

Figure  4-9  presents  a  dielectric  loaded  waveguide  problem.  We  have  solved  this  problem  and 
compared  our  results  with  the  approximate  solutions  of  Marcuvitz  ,  who  estimated  the  error  in 
his  approximation  to  be  less  than  a  few  percent.  This  comparison  is  presented  in  Table  4-2. 


116 


Magnitude  and  phase  of  reflection  coefficient  of 
H-plane  step  discontinuity 


Magnitude  and  phase  of  reflection  coefficient  of 
E-plane  step  discontinuity 


reflection  coefficient  of 
[-plane  discontinuity) 
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Figure  4-7:  Single  half-round  inductive  obstacle 
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Figure  4-8:  Waveguide  bending 


Table  4-1:  Tbe  reflection  and  transmission  coefficients  R  and  T,  relative  to 

the  reference  planes  t,  were  calculated  by  present  method  for 

H-plane  corners  with  9  =  30°,  9  =  60°,  9  =  90° 

and  compared  with  experimental  data  from  Marcuvitr24. 


X 

1  - 

90 

#  - 

80 

I  - 

00 

a 

experimental 

present 

experimental 

present 

experimental 

present 

data 

method 

data 

method 

data 

method 
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R—0.020.j0.022 

R-0.130-fjO.057 

R— 0.144+j0.054 

R— 0.053+ j0.4 17 

R— 0.041+j0.430 

0.762 
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T— 0.000+ j0.114 

T— 0.808+ j0.084 
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R—0.02S-j0.031 
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R— 0.236+ jO.202 

0.6724 
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Figure  4-8:  Dielectric  post  waveguide 


Table  4-2: 


Dielectric  post  in  rectangular  guide  (H]0  mode)  our  results 
compared  with  the  approximate  analytic  solutions 
calculated  from  the  formula  in  Marcuvit*24 
for  different  dielectric  constants.  (X/a  **  1.4) 


Approximate 

Analytic 

solution 

Present  method 

‘r 

mesh  1 

matrix  size  146 

mesh  2 

matrix  size  260 

4.0 

R  =  -0.025- jO.  155 

T  =  0.074  -  jO.157 

R  = -0.021  -  jO.  142 

T  =  0.982 -jO.  130 

R  =  -0.022  -  j0.148 

T  =  0.077  -  jO.  140 

0.0 

R  =  -0.101  -  j0.381 

T  *  0.809  -  j0.396 

R  =  -0.172-  jO.379 

T  «  0.832  -  j0.371 

R  ■>  -0.175  -  j0.378 

T  —  0.824  -  j0.382 

16.0 

R  =  -0.588  -  jO.  490 

T  *  0.431  -  jO.501 

R  *  -0.528  -  jO.SOO 

T  *  0.474  -  j0.498 

R  =  -0.534  -  j0.404 

T  -=  0.466  -  j0.503 

a 

R  *  -0.882  -  j0.313 

T  »=  0.117  -  j0.331 

R  =  -0.851  -  j0.354 

T  *  0.150  -  j0.360 

R  =  -0.850  -  j0.339 

T  =  0.140- j0.356 
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SECTION  5  THE  TWO-BOUNDARY  APPROACH  FOR  OPEN 
SCATTERING  PROBLEMS 

6.1.  INTRODUCTION 

An  important  practical  application  of  electromagnetic  field  theory  is  the  prediction  of  the  radar 
cross-sections  of  bodies  that  lie  in  unbounded  space  and  are  illuminated  by  an  electromagnetic 
wave.  Radar  cross-sections  provide  a  measure  of  the  amount  of  energy  scattered  by  a  body  in  a 
given  direction  and  may  be  calculated  if  the  electric  field  around  the  body  is  known. 

In  the  waveguide  junction  problem  of  Section  4,  the  governing  equation  is  the  Helmholtz 
equation  and  the  geometry  extends  to  infinity  only  one  dimension.  Using  the  same  governing 
equation  but  extending  to  infinity  in  two  dimensions,  open  cylindrical  scattering  problems  can  be 
treated  as  a  modification  of  waveguide  junction  problems.  The  discontinuities  in  the  waveguide 
which  cause  the  perturbation  of  the  fields  is  replaced  by  scatterers  in  free  space  that  deflect  the 
incident  wave.  Further,  the  reflected  and  transmitted  waves  we  computed  in  Section  4  are  now 
the  scattered  fields  from  the  structure. 

The  cross  section  of  a  cylinder  with  refraction  index  q  is  shown  in  Figure  5-1.  Define  the 
incident  field  E"  and  H1  as  the  field  with  the  cylinder  absent  and  the  scattered  field  E®  and  Hs  as 
the  difference  between  the  total  field  with  the  cylinder  present  (E1,^)  and  the  incident  field,  that 

is, 

£*  =  £*-£*  H*  =  H*  -  H* 

By  assuming  that  the  cylinder  is  infinitely  long  and  that  the  incident  wave  propagates  in  the 
plane  of  the  figure,  the  scattering  problem  is  reduced  to  two  dimensions.  When  the  length  of  the 
cylinder  is  very  large  compared  to  the  wavelength,  this  assumption  will  not  contribute  too  much 
error  in  the  analysis.  However  if  the  length  of  the  cylinder  is  of  the  same  order  as  the  wavelength 
then  a  three  dimensional  analysis  is  needed . 
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Figure  6-1:  The  cross  section  of  a  scatterer 


In  this  Section,  we  treat  only  the  case  of  a  TE  plane  incident  wave  having  the  electric  field  in 
the  i  direction.  For  the  other  cases,  such  as  incident  TM  waves  or  cylindrical  wave,  the  analysis 
follows  in  much  the  same  way.  From  Maxwell  equations,  the  total  electric  Held  with  the  cylinder 
present  must  satisfy  the  Helmholtz  equation  in  all  space 

(V2  +  it02  »f)  £*  =  0 

kQ  is  the  wavenumber  in  free  space  and  tj  is  the  refraction  index  of  the  material.  Furthermore,  the 
scattered  field  will  die  out  as  we  go  away  from  the  cylinder,  the  boundary  conditions  at  infinity  is 
El  ==  E1.  Finally  we  will  adopt  the  polar  coordinates  (r-0)  system  to  analyze  open  scattering 
problems  since  here  two  dimensions  extend  to  infinity. 

6.2.  FORMULATION  AND  APPLICATION  OF  TWO-BOUNDARY  APPROACH 

6.2.1.  MODAL  ANALYSIS 

For  simplicity,  assume  that  the  incident  wave  is  a  TE  plane  wave  that  propagates  in  the  x 
direction  and  that  the  geometry  of  the  scatterer  is  symmetric  with  respect  to  the  x  axis.  We 
define  a  circle  r j  with  radius  Tj  to  serve  as  the  boundary  of  the  finite  element  solution  region, 
and  call  it  the  "picture  frame  circle".  A  second  circle  /"2  with  smaller  radius  r2  is  defined  to  be 
the  mode  matching  boundary  and  is  called  the  "matching  circle",  as  shown  in  Figure  5-2. 
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Figure  5-2:  The  picture  frame  circle  and  the  matching  circle 

The  reasons  for  choosing  circles  as  the  picture  frame  and  matching  boundaries  are: 

1.  The  circle  is  the  simplest  shape  to  be  used  as  a  boundary  in  polar  coordinates. 

2.  The  analytical  solution  for  the  exterior  region  is  given  by  Hankel  functions  in  the  r 
direction.  Computing  the  Hankel  functions  is  an  expensive  task  compared  to  other 
numerical  procedures.  By  using  circles  for  the  picture  frame  and  matching  boundaries 
we  only  need  to  compute  the  Hankel  function  once  at  each  radius. 

3.  Since  Hankel  functions  oscillate  and  change  many  orders  of  magnitude,  the  numerical 
solutions  would  be  unstable  if  the  boundaries  involved  radii  of  different  magnitudes. 

For  a  circular  boundary,  only  one  radius  is  involved  and  this  leads  to  stable  numerical 
procedures. 

The  analytical  solution  for  the  scattered  Held  in  the  region  r  >  r2  can  be  found  by  using  the 
method  of  separation  of  variables1.  The  result  is 

nt{r,0)  =  Hfar)  coe{i9)  (5.1) 

due  to  the  symmetry,  only  the  cos  terms  are  introduced  in  (5.1).  kQ  is  the  wavenumber  in  the 
free  space,  and  H.  is  the  Hankel  function  of  second  kind  of  the  order  i.  If  we  use  (M+l)  modes  to 
approximate  the  scattered  electric  Held,  then,  together  with  the  incident  electric  field,  the 
analytical  expression  for  the  total  field  in  the  region  r  >  r2  is 

M 

E?(r,6)  =  exp{-jkoreo«0)  +  ^  ap,lr<e)  r  -  r2  (5-2) 

i=0 

where  expf-jkfotd)  is  the  incident  electric  field  and  a.’s  are  the  coefficients  that  need  to  be 
determined.  Note  that,  in  general,  the  coefficients  a_  are  complex  since  all  the  functions  are 
complex  valued  in  (5.1). 
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Since  equation  (5.2)  holds  for  r  >  r 2,  it  gives  the  electric  field  at  the  boundary  Fj  as 


M 


Ef(rvB)  =  ex^-^fjcostf)  +  £  a,'»,(rr^ 


i-0 


We  may  approximate  the  term  cxpf-jkj-fosd)  by  its  (M+l)  Fourier  components  as 


M 


where 


exp{-jk0rxcoeB)  —  cQ  +  £  exos(»«) 
i— 1 


J*'  txp(-]\rxco*(B))d6 
l  /-a* 

e,  =  “  /  exp(-y*0r1co«(tf))*co«(itf)  dtf 
Thus,  equation  (5.3)  can  be  written  as 

E*(ryB)  =  {«0H0(*0ri) +  co}  +  £  {  +  c.)  *  coa(,’fl) 


(5.3) 


(5.4) 


From  equation  (5.4),  we  set  up  (M+l)  boundary  value  problems  for  the  field  patterns  X.(r,0) 
inside  the  finite  element  solution  region.  X.(r,0)  is  obtained  by  using  the  finite  element  method  to 
solve  the  Helmholtz  equation  with  the  boundary  conditions  X.(r1(0)  =  cos(itf)  on  the  boundary  Fj 
within  the  finite  element  solution  region.  Once  the  (M+l)  field  patterns  are  known,  by  using  the 
superposition  theorem,  the  total  field  in  the  finite  element  solution  region  can  be  expressed  as 
follows 

TZ  ‘  (r,tf)  =  {a0H0(k0r,)  +  cQJ  \(r,B)  +  £  {a .  HfarJ  +  cj  X^r.fl)  (5.5) 

i=l 

Again,  both  equations  (5.2)  and  (5.5)  are  valid  in  the  region  r2  <  r  <  r^,  so  that  we  have  two 
solutions  for  the  field  on  the  matching  circle  F2 

M 

E?(r2,B)  mm  cxp(-jk0r2eoeB)  +  £  a.  ifa.B)  (5.6) 


i-0 


and 


(5.7) 


i— 0 


The  last  step  is  to  match  these  two  solutions  on  the  circle  F2,  from  which  the  complex 
coefficients  a’s  can  be  solved. 

i 
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5.2.2.  MATCHING  FORMULATION 

Matching  the  two  solutions  is  performed  by  minimizing  the  square  of  the  difference  between 
(5.6)  and  (5.7)  on  the  matching  circle.  From  (5.6)  and  (5.7)  we  construct  an  L2  norm  of  the 
difference  between  the  numerical  solution  and  the  analytical  solution  on  the  matching  circle 


L2  =  f0*  \T:\rvo)-Ef(r2,«)\2de 

and 

J:  \r2,9)  -  E*{rv0)  |  2  =  {/C(*)  +  £  a .  Q  .(*)}•{*>)  +  £  a.  Q*(<?)} 

i=0 


i=0 


where 


M 


K(9)  =  exp{-jk0r2co8$)  -  ^  c,  \(r2’*) 


i=0 


QJM  =  H(kQr2)co8(i9)  -  Ht(kori)\t(r2,9) 

where  the  superscript  represents  the  complex  conjugate.  By  writing  aj  =  x.  +  j  y;  and 
minimizing  L22  with  respect  to  x.’s  and  y.’s,  the  following  equations  are  obtained 

E  «>+E  ©>  -  -/**<*•  «,} ■» 


1=0 


M  M  _2jf  /*2jf 

£*,/,  Q')d9+Ylyif<i  R'{QiQ'i)de-  J0  Irn{K  Qt)d9 

«=0  »=0  J 


1=0 


M 


1  =  0....M 


(5.8) 


Finally  from  (5.8)  we  can  solve  for  x.’s  and  y.’s  and  obtain  the  coefficients  a.’s. 


5.2.3.  FAR  FIELD  PATTERN 


Once  the  complex  coefficients  a’s  are  known,  the  scattered  field  in  the  exterior  region  is 
obtained  as 


M 

£*M)  =  £  MiM) 

1=0 


£  r2 


(5.9) 


At  large  distances  (r  >  >  wavelength  ),  the  scattered  field  in  terms  of  (5.9)  can  be  written  as 
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2/  N 
E?  =  (^-H1/2  exp(-jkQr)»{a0  +J^3*  ai  co*(«'*)} 


It  is  customary  to  define  the  far  field  scattering  pattern  by  the  equation 
N 

g(ff)  =  aQ  +  £  J*  a.  coe(iS) 
i— l 


6.3.  NUMERICAL  RESULTS 

Shown  in  Figure  5-3  is  the  geometry  of  two  parallel  circular  cylinders  with  radii  0.2  wavelength 
and  0.1  wavelength  separated  by  0.4  wavelength.  The  refraction  indexes  are  2.0  for  both 
cylinders.  This  problem  has  been  analyzed  by  Mei11.  The  radius  of  the  picture  frame  circle  is  0.6 
wavelength  as  also  shown  in  the  figure.  Figure  5-4  shows  two  different  meshes  used  to  solve  this 
problem:  Figure  5-4  (a)  contains  165  unknowns  and  Figure  5-4  (b)  contains  287  unknowns. 

6.3.1.  INCIDENT  ANGLE  fllnc  =  0* 

When  a  plane  wave  is  incident  with  angle  0.bc  =  0°  from  the  x  axis,  the  solution  should  be 
symmetric  with  respect  to  the  x  axis.  Figure  5-5  presents  a  plot  of  lines  of  equal  values  of  Ez  at 
time  phase  0°,  calculated  by  using  N  =  4  and  r2  =  0.4  wavelength  in  the  program.  The  jagged 
lines  in  this  figure  are  a  result  of  imperfections  in  the  plotter  and  do  not  represent  discontinuities 
in  the  solutions  which  were  smooth.  The  resultant  amplitude  of  the  far  field  pattern  from  the 
program  compared  with  solutions  obtained  by  Mei11  are  shown  in  Figure  5-6  for  the  two  meshes 
(a)  and  (b). 

6.3.2.  INCIDENT  ANGLE  0|ne  =  80* 

Shown  in  Figure  5-7  are  the  plots  of  equal  values  of  Ej(  for  an  incident  angle  of  9-me  —  90  °  . 
The  parameters  used  in  the  program  are  :  N  *  4,  r2  0.4  wavelength,  and  the  mesh  used  is  the 
same  as  Figure  5-4(b).  Figure  5-8  is  a  comparison  of  the  amplitude  and  phase  of  the  far  field 
pattern  obtained  by  the  present  method  and  that  obtained  by  Mei11. 


Figure  5-3:  Two  parallel  circular  cylinders  with  radii 
0.2  wavelength  and  0.1  wavelength 
separated  by  0.4  wavelength 


(b) 


Figure  6-4:  Finite  element  meshes  for  the  problem  in  Figure  5-3 

(a)  165  unknowns  (b)  287  unknows 
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(b) 


Figure  5-5:  Plots  of  equal  values  of  with  time  phase  0° 

(a)  for  the  mesh  in  Figure  5-4  (a) 

(b)  for  the  mesh  in  Figure  5-4  (b) 
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ittern  for  6.  =0 


Figure  5-0:  Amplitude  of  the  far 


isnrr 


(a)  time  phase  ==  0  * 

(b)  time  phase  =  90  * 


Figure  6-8:  Far  Held  pattern  (a)  amplitude  (b)  phase 
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SECTION  6  THE  TRANSFINITE  ELEMENT  METHOD  FOR 
ELECTROSTATIC  PROBLEMS 


Described  in  this  chapter  and  in  Section  7  is  a  novel  coupling  scheme  for  the  open  problems 
called  the  "transfinite  element  method*.  The  purpose  of  this  name  is  to  denote  the  fact  that 
analytical  basis  functions  are  employed  in  the  infinite  exterior  region.  In  this  chapter,  we  take 
the  electrostatics  as  a  simple  case  to  develop  the  theory.  In  Section  7,  we  return  to  the  open 
scattering  problems  but  using  the  transfinite  element  method. 

tt.l.  GENERAL  DESCRIPTION  OF  THE  TRANSFINITE  ELEMENTS  METHOD 

The  solution  procedure  for  applying  the  transfinite  elements  method  is 

1.  Derive  a  functional  either  by  using  the  Rayleigh-Ritz  method  or  by  using  Galerkin's 
method  for  the  problem  being  solved.  The  functional  derived  will  be  an  integral  over 
the  entire  problem  region. 

2.  Create  a  picture  frame  that  separates  the  problem  region  into  two  parts  /?.  and  Q 
The  region  interior  to  the  picture  frame  is  (2.  and  the  region  exterior  to  the  picture 
frame  is  Let  r  be  the  boundary  between  them.  The  finite  element  method  is  used 
in  the  interior  region,  so  the  picture  frame  must  enclose  all  inhomogeneities  and  other 
quantities  that  must  be  modeled  by  the  finite  element  method.  The  functional  integral 
over  the  entire  region  is  thus  separated  into  two  parts:  one  integral  over  the  interior 
region  17.  and  the  other  over  the  exterior  region  f?e- 

3.  Because  the  exterior  region  is  homogeneous,  we  can  create  a  complete  function 
space  Af  for  all  of  the  functions  that  satisfy  the  boundary  value  problem  in  the 
exterior  region  by  using  the  superposition  theorem  and  the  method  of  separation  of 
variables.  The  basis  functions  in  thus  found  will  be  exact  solution  of  the  problem  in 
the  exterior  region.  It  follows  that  the  integral  of  the  functional  over  17 t  is  tero  and 
that  only  the  integral  over  the  interior  region  Q .  must  be  computed. 

4.  Apply  the  standard  finite  element  method  in  the  interior  region.  The  basis  functions 
employed  can  be  the  conventional  Lagrange  interpolation  polynomials.  The  function 
space  A  constructed  from  these  polynomials  represent  the  the  basis  vectors  for  the 
function  space  for  the  numerical  solution  in  the  interior  region. 
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5.  Impose  continuity  conditions  of  the  solution  between  the  function  space  A  and  the 
function  space  Af  along  the  boundary  r.  This  may  be  accomplished  simply  by  setting 
the  values  from  both  solutions  equal  at  the  boundary  nodes.  Note  that  usually,  the 
number  of  unknowns  used  to  approximate  the  solution  in  the  exterior  region  is  less 
than  the  number  of  nodes  on  the  boundary  F.  The  final  function  space  A  for  the 
numerical  solutions  in  the  whole  region  will  be  the  union  of  the  function  space  A  and 
A  restricted  by  the  continuity  conditions  imposed  on  the  boundary  F. 

6.  Finally,  apply  the  stationarity  property  of  the  functional  on  the  function  space  A  to 
obtain  a  matrix  equation  to  be  solved  for  the  coefficients  in  the  approximation. 


•.2.  FORMULATION  OF  TRANSFINITE  ELEMENTS  IN  THE  ELECTROSTATIC 
PROBLEMS 

6.2.1.  FUNCTIONAL 

The  particular  problem  with  which  we  will  be  concerned  here  is  the  solution  of  the  Poisson 
equation 

—  V  •  t  V  ip  —  a  in  IP 

r  *  *  oo 

subject  to  the  boundary  condition 

,  „  (6.1) 
\p  —  0  at  oo 


The  problem  region  12  for  (6.1)  is  of  course  unbounded. 

According  to  the  Rayleigh-Ritz  method,  the  solution  of  (6.1)  may  be  obtained  by  minimizing  the 
following  functional 


—  f  trV  rp  ip  dlP  —  2  f  g  ip  dIP 

J  00  00 


(6.2) 


By  using  a  circle  F  with  radius  Tj  to  enclose  all  inhomogenities  and  sources  ,  we  separate  l?o 
into  two  regions:  an  interior  region  12.  and  the  exterior  region  J?(,  as  shown  in  Figure  6-1. 


The  integral  in  (6.2)  can  be  integrated  separately  over  12.  and  12 1  as  follows 

V  ip dI7  —  2  ^  grp dI2  +  V  \p*^  \p dTP 


(6.3) 


By  Green’s  theorem  we  may  rewrite  this  as 
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Figure  6-1:  The  separation  of  the  problem  region 

interior  region  J7.  exterior  region  and  the  boundary  f 


dF~  (6.4) 

The  reason  for  converting  the  functional  from  (6.3)  to  (6.4)  will  become  clear  in  the  following 
section. 


F\i’)  =  Jn(r^  ^ dO -  2 grp dn - 


6.2.2.  THE  FUNCTION  SPACE 


Note  that  the  exterior  region  i?t  now  is  homogeneous  and  source  free.  If  the  net  charge 
contained  in  the  interior  region  /?.  is  zero,  the  function  space  Af  obtained  by  using  the  method  of 
separation  of  variables  is 


=  ~i  *  (  Voa(’*)  +  V,n(‘*) )  r  —  rl  } 

i—l  r 


(6.5) 


where  M  is  the  number  of  terms  used  to  approximate  the  solution  in  the  region  f?t,  and  a  and  b; 
are  scalar  numbers.  If  the  net  charge  in  f?.  is  not  zero,  a  logarithmic  term  that  vanishes  at  a 
certain  reference  point  must  be  included  in  the  expansion. 


We  subdivide  the  interior  region  f2.  into  second  order  triangular  finite  elements  and  use 
Lagrange  polynomials  as  the  interpolation  functions.  The  function  space  A  for  the  solution  in  the 
n.  is 

i 
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I  +  of  r<,,}  (6.6) 

6—1  1—1 

where  we  have  separated  the  interpolation  polynomials  into  o.1  for  the  nodes  in  ft  but  not  on  the 

r  i  r 

boundary  F  and  a.‘  for  the  nodes  on  the  boundary  r.  Co rre ponding  to  a/  and  o  '  are  the 

»  r 

solution  values  tP.  and  tfy  ,  respectively.  This  subdivision  will  help  in  imposing  continuity 

conditions  in  the  procedure. 

The  combined  function  space  A  m  A^  U  A*  as  discussed  in  section  6.1  form  the  new  function 
space  A 

A- {v-l  l*  m  n*  ,  i>  6  C?)  (6.7) 

1  IV'  in  ft  > 

Continuity  of  the  field  is  imposed  by  setting 
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where  (rlt8.)  is  the  polar  coordinates  for  the  boundary  point  V>.r,  and 
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6.2.3.  THE  BOUNDARY  INTEGRAL 

Since  the  functions  in  A  are  analytical  solutions  of  the  Laplace  equation,  the  integral  over  Q  in 
the  functional  (6.4)  is  sero  for  any  choice  of  the  coefficients  a.  and  b..  Thus  the  functional  reduces 
to  the  following  form 


A  I  ( 

=  [  t  vtf  •  <m -  2  [  gtfdn-  f  t'—dr 

Jni  Jnt  Jr  dr 


(6.8) 


Thus  by  a  judicious  choice  for  the  function  space  A,  the  functional  F(0)  is  reduced  from  integrals 
over  the  entire  unbounded  region  into  integrals  over  the  interior  region  plus  a  boundary  integral 
term.  This  reduction  has  also  be  exploited  by  Silvester  and  Hsieh9,  although  their  procedure  for 
modeling  the  exterior  region  was  quite  different. 

The  boundary  integral  term 


r 

{/- 


dr 


in  (6.8)  is  evaluated  simply  by  analytical  integration.  From  (6.5),  we  have 
H *  M 

i—i 


Subsitituting  (6.10)  and  (6.5)  into  B  gives  the  result 
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(6.11) 


where 
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6.2.4.  THE  RESULTING  MATRIX  EQUATION 


From  equation  (6.8),  the  nodal  values  V1'  can  be  written  as  a  connection  matrix  times  a  column 
rector 


where  NIN  is  an  N  by  N  identity  matrix,  and 


0  8  R 


Substituting  (6.11)  and  (6.12)  into  the  functional  (6.9)  results  in 
F[x)  =  x^  x  —  2  C  i  +  D  A 

where 


(6  13) 


where  the 


below  a  variable  indicates  a  row  vector,  and  T  represents  the  transpose  of  a  matrix. 
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Finally  minimizing  the  functional  F(x)  with  respect  to  the  variables  x  gives  the  resulting  matrix 
equation 


sil  sir  0R 

ReTsrr0R 

R0TsrI  +  d 


g  I 


(6  1-0 

Notice  that  the  matrix  on  the  left-hand-side  is  real,  symmetric  and  positive  definite  and  can 
therefore  be  solved  very  efficiently  using  the  incomplete  Choleski  pre-conditioned  conjugate 
Gradient  Method26.  Also  note  that  if  the  number  of  coefficients  for  the  exterior  region  is  less  than 
the  number  of  nodes  on  the  boundary  f,  as  is  often  the  case,  that  the  total  number  or  unknowns 
is  reduced  by  applying  the  transfinite  element  procedure. 

#.3.  NUMERICAL  RESULTS 


6.3.1.  PARALLEL-PLATE  CAPACITOR 


The  problem  of  a  square,  parallel-plate  capacitor  was  chosen  as  an  example  to  illustrate  the 
method.  A  cross-section  of  the  geometry  is  shown  in  Figure  6-2  and  represents  two  infinitely  long, 
thin  conducting  plates.  This  problem  is  of  particular  interest  since  it  contains  a  singularity  at  the 
edge  of  the  capacitor  plate.  The  mesh  used  for  this  problem  is  shown  in  Figure  6-3  and  a  plotof 
the  equal  potential  lines  is  shown  in  Figure  6-4.  As  is  evident  from  the  plot,  the  behaviors  of  the 
solution  inside  the  picture  frame  is  the  same  as  one  would  expect  if  all  space  were  modeled 

6.3.2.  CIRCULAR  CONDUCTORS 


The  second  example  is  two  infinitely  long  circular  conductors  both  with  radius  A  and  separated 
by  distance  D,  the  ratio  A/D  is  equal  to  0.25.  The  exact  value  of  the  capacitance  per  unit  length 
of  this  system  is  given  by  the  formula27 


eo«h  1{D/2A) 
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la  this  case  the  value  of  C  =  2.3855  <Q.  Due  to  the  symmetry  of  the  problem,  we  used  only  one 
half  of  the  problem  region.  The  geometry  of  the  problem  and  the  semi-circle  which  served  as  the 
boundary  F  is  shown  in  Figure  5-5;  the  mesh  and  the  resulting  equi-potential  line  plot  are  in 
Figure  6-6  and  Figure  6-7,  respectively.  The  capacitance  we  obtained  is  2.3783  <0>  an  error  of 
0.39c. 


6.3.3.  FOUR  LAYERS  CAPACITOR 

The  last  electrostatics  example  is  two  infinitely  long  rectangular  cylinders  filled  with  four 
different  dielectrics.  The  relative  dielectric  constants  are  2.0,  3.0,  4.0,  5.0  for  layer  I,  II,  III,  IV;  as 
shown  in  Figure  6-8.  A  simple  approximation  which  treat  these  four  capacitors  in  series  gives  the 
C  «=  3.97351  <Q  per  unit  length.  Analysis  by  the  transfinite  element  method  gives  a  value  5.044 
<Q.  Based  on  the  appearance  of  potential  plot  in  Figure  6-10,  we  expect  that  the  transfinite 
element  analysis  provides  the  better  result. 


Figure  8-4:  The  equal  potential  lines  plot.  The  potential  on 
the  top  conductor  is  +1  volt  and  the  potential  on  the  bottom  conductor  -1  volt 


igure  ft-5:  Circular  cylinders  with  D 


Figure  Mesh  for  example  2 
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Figure  6-7:  The  equal  potential  lines  plot  for  example  2.  The 
potential  on  the  circular  conductor  is  1  volt  and  the  potential  on  the 
ground  plane  is  0  volt 


Figure  8-8:  The  geometry  of  Tour  layers  capacitor, 

with  relative  dielectric  constants  2.0,  3.0,  4.0,  5.0 


Figure  8-0:  The  mesh  for  the  four  layers  capacitor 
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SECTION  7  THE  TRANSFINITE  ELEMENTS  METHOD  FOR  OPEN 
SCATTERING  PROBLEMS 


7.1.  FORMULATION  IN  TERMS  OF  THE  TRANSFINITE  ELEMENTS  METHOD 

7.1.1.  WEAK  FORMULATION  OF  GALERKIN’S  METHOD 

As  mentioned  in  Section  2,  variational  principles  can  be  obtained  either  by  using  Galerkin’s 
method  or  by  using  the  Rayleigh-Ritz  method.  It  worth  pointing  out  that  Galerkin’s  method  is 
more  general  than  Rayleigh-Ritz  method  since  it  allows  non-self  adjoint  operators  to  be  modelled. 

By  applying  Galerkin’s  method  as  discussed  in  Section  2,  the  following  bilinear  form  is 
obtained 

B(u,u)=  f  u*  (V2u  +  kQ2tf  u)  dn  (7.1) 

J  00 

where  1?^  is  the  problem  region  and  represents  the  complex  conjugate.  As  shown  in  Figure  6-1, 
the  circle  /"with  radius  ^  separates  the  problem  into  i7  and  f?e.  From  the  Green’s  theorem,  the 
bilinear  form  B(v,u)  in  (7.1)  can  be  written  as 

d%L 

(Vv  *Vu  —  v  kfau)  dn  —  j v* ~  dr  —  V  (V2u  +  fr2ti)  dn  (7.2) 


7.1.2.  THE  FUNCTION  SPACES 

The  trial  function  space  for  u  and  the  testing  function  space  for  v  are  set  equal  in  the  analysis 
here.  The  scattered  field  E*  in  the  region  i?e  can  be  approximated  by 

+ YL  H,{korHafO9{i0) + (7-3) 

•=i 

Combining  this  with  the  incident  electric  field  ezpf-jkfcoeO)  gives  the  function  space  Af 
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Figure  7-1:  Interior  region  Q.,  Exterior  region  l?e 
and  the  boundary  r 

\  =  {’/''I  tl>e^cxp(-]TcQrco80)+aoHo{kor)+J2  HtikQr)* {afosW+bfintfe) } }  (7-4) 

i—1 

The  function  space  A  in  the  interior  region  1?.  is 

*■-£;  ♦/«/+£  *[«[ 

i«=l  i»=l 

Where  the  same  notation  is  used  here  as  in  Chapter  6.  Continuity  of  fields  at  the  boundary  r 
results  in  the  conditions 
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(7.5) 

Note  that  the  matrix  H  ia  diagonal.  The  function  space  A  is  the  union  of  Af  and  A.  restricted 
with  the  continuity  conditions  (7.5). 

Since  u  and  v  are  both  in  the  function  space  A  they  must  satisfy  the  Helmholti  equation  in  the 
exterior  region  exactly,  then  the  bilinear  form  B(v,v)  simplifies  to 

EKv.u)  -  Ja  (Vu*.Vu  -  vfy  u)  dn-  j/^dr  (7.6) 

The  functions  u  and  v  can  be  expressed  as 
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N  P 

»— i  i-i 


£*nc  +  aQ»HQ(kQr)  +  £  Hn(k0r)*{  ancot(ne )  +  *„«'"("*)) 


£aM+I>^  ,n/?. 

"“{'"I  *'*1  M  (77) 

£*’ne  +  c0«//0(i:0r)  +  £  Hn(k0r)»(  eneoe{n«)  +  d«n(n*))  in  f?e 


7.1.3.  BOUNDARY  INTEGRAL  TERM 


Note  that  since  Fis  a  circle  with  radius  r^  the  boundary  integral  term  simplifies  to 


r  *ou  fit  .  du 

lrvTndr~rfj<,  ’  " 


From  equation  (7.7),  we  have  that 


«*(»v#) — E  ‘ne(«v*)+‘X(Vi)+E  //lVi)'(c«,CM*n>*+<,»"n*nJ*) 


3u  a£*nc  M  f  \ 

jfirV9) - ar^rrtf)+VWl)+2I  \(Vl)#(anco,(ntf)+6»*,n(nfl)) 


where  hlk.r.)  =  — - — |  .  Approximating  E*ne  and  — — —  by  (2M+1)  Fourier  components 

•  u  *  dr  Vi  dr 


£*ne(»V*)  —  e0+  {  tnC09(n9)  +  /B«'n(n^)  J 

»■»! 

a£**e  w  , 

~  ff0+  £  \  Jnco«(n«)  +  ln»in(nO) } 

Thus  substituting  (7.9)  into  (7.8)  ,  integrating,  and  using  the  orthogonality  property  of 
trignometric  functions,  the  boundary  integral  (7.8)  assumes  the  matrix  form 
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7.1.4.  A  MATRIX  EQUATION  FOR  SCATTERING 

From  equations  (7.5)  and  (7.10),  we  see  that  the  functional  B(v,u)  can  be  expressed  in  the 
matrix  form 


where 


(7.11) 
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Vo  •Vo>  -  a.  kfr  a. )  dn 


To  obtain  a  weak  solution  for  u  ■*  4/,  B(v,4>)  —  0  for  V  v  t  A.  This  condition  results  in  the 
following  complex  matrix  equation  to  be  solved  for  the  coefficients  7  '  and  X 
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(7.12) 

Equation  (7.12)  can  be  simplified  to  give  the  final  results 
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7.2.  NUMERICAL  RESULTS 

The  validity  and  versatility  of  the  techniques  presented  here  are  demonstrated  by  the  following 
examplpes. 

7.2.1.  A  CIRCULAR  CYLINDER  WITH  ITS  CENTER  OFFSET  FROM  THE  ORIGIN 

Numerical  solutions  are  generated  by  the  transfinite  element  method  with  the  center  of  the 
boundary  circle  F  at  the  origin,  while  the  center  of  the  scatterer,  which  has  radius  0.3  X,  is  offset 
by  0.2  X  from  the  origin,  as  shown  in  Figure  7-2.  Figure  7-3  is  the  mesh  used  in  the  program, 
which  corresponds  to  187  unknowns.  Placing  it  off-center  in  this  analysis  provides  a  good  test  of 
the  validity  of  the  method. 

I.  Figure  7-4  is  the  far  field  pattern  |g(^)|  for  the  refraction  index  of  the  scatterer  equals  to  2.0, 
excellent  agreement  between  numerical  and  exact  solution  exists.  The  exact  solutions  in  Figure 
7-4  are  obtained  by  using  the  method  of  separation  of  variables11.  The  equal  electric  field  plots 
for  different  time  phases  are  shown  in  Figure  7-5.  The  elements  outside  the  circle  in  Figure  7-3 
are  used  to  interpolate  the  field  values  in  the  exterior  region.  The  electric  field  outside  the 
boundary  circle  Fhas  been  evaluated  by  using  equation  (7.4). 

II.  Shown  in  Figure  7-6  are  the  intensity  plots  and  the  far  field  patterns  |g(  v'’)|  for  refraction 
indeces  of  2.0+j2.0,  2.0+j20.0  and  2.0+j200.0.  The  strange  curves  within  the  shaded  area  are  due 
to  the  fact  that  the  field  intensities  are  almost  zero  there.  Unfortunately,  analytic  solution  for 
these  lossy  cases  do  not  exist  for  comparison. 

7.2.2.  PLANE  WAVE  SCATTERING  BY  A  ELLIPTIC  CYLINDER 

Figure  7-8  presents  the  far  field  pattern  |g(V>)[  for  the  elliptic  cylinder  scatterer  in  Figure  7-7 
with  a  refraction  index  of  2.0.  The  results  calculated  here  are  compared  with  solutions  from 
Marin7  and  again  varify  the  accuracy  of  the  transfinite  element  approach.  Propagation  patterns 
at  different  time  phases  are  shown  in  Figure  7-9. 

7.2.3.  SCATTERER  WITH  NONCONVEX  NONSYMMETRIC  CROSS-SECTION 

The  cross-section  of  a  kidney-shaped  scatterer  is  shown  in  Figure  7-10(a)  where  the  dimension  of 
the  maxmimum  object  radius  is  0.32  X.  The  incident  radiation  is  a  plane  wave  with  wavelength  1 
cm,  and  the  direction  of  the  incident  wave  is  45°  from  the  positive  x-axis.  Three  different 
materials  with  refraction  indexes  7,  6+j  200  and  j  20000  are  investigated,  which  are  correspond 
approximately  to  glass,  salt,  and  carbon.38'20 
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Figure  7-10(b)  provides  the  finite  element  mesh  used  in  the  program. The  results  of  these 
calculations  are  shown  in  Figure  7-11  and  are  arranged  in  order  top  to  bottom  of  increasing 
conductivity.  The  graphs  on  the  left  side  of  Figure  7-11  are  plots  of  the  intensity  of  the  total 
wave.  The  graphs  on  the  right  side  of  Figure  7-11  display  |g(0)|  versus  ^  for  these  three  materials 
along  with  a  comparison  with  the  results  of  Marin7. 


Figure  7*3:  Mesh  size  187  unknowns 


(a)  time  phase 


(d)  time  phase  00 


ure  7-6t  Equal  line  plots  of  E  field  for  different  tune  phases 
(a)  time  phase  0  (b)  time  phase  30 

(c)  time  phase  00  (d)  time  phase  90  " 


(c)f)  —  S.O  +  jSOO.O 


Figure  7-8:  The  intensity  plots  and  the  far  Held  pattern  plots  for  various  refraction  indexes 

(a)  2  +  j  2.  (b)  2  +  j  20.  (c)  2  +  j  200. 
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(b)  Matrix  size  214 


Figure  7-7:  Dielectric  elliptic  cylinder  and  the  corresponding  mesh 
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(b)  mesh  -  Matrix  size  314 


Figure  7-10:  The  geometry  (a)  and  mesh  (b)  for  a  kidney-shaped 
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APPENDIX  F  -  AN  ADAPTIVE  SPECTRAL  RESPONSE  MODELING 
PROCEDURE  FOR  ELECTROMAGNETIC  SCATTERING 


ABSTRACT 

An  adaptive  scheme  is  proposed  to  generate  the  spectral  response  of  waveguide  junctions  in  minimum 
computation  time.  The  procedure  uses  the  newly  developed  transfinite  element  method  to  determine  the 
fields  in  junctions  at  a  few  adaptively  selected  frequencies  and  then  employs  these  solutions  to  generate 
the  spectral  response  throughout  the  frequency  range  of  interest.  In  typical  problems,  the  method  con¬ 
verges  in  five  or  six  iterations  to  the  full  spectral  response  evaluated  at  one  hundred  points.  We  show  by 
solving  example  problems  that  the  new  procedure  is  orders  of  magnitude  faster  than  the  alternatives. 

INTRODUCTION 

Microwave  circuits  in  use  today  often  employ  planar  geometries  that  may  be  represented  mathemati¬ 
cally  as  an  Alport  microwave  junction.  An  N-port  microwave  junction  is,  in  general,  defined  as  a  struc¬ 
ture  that  consists  of  an  arbitrarily  shaped  cavity,  with  or  without  dielectrics,  and  has  N  rectangular  ports 
coupling  in  and  out  of  the  cavity.  A  number  of  studies  have  been  made  of  microwave  junction 
problems1,2’3,4.  However,  the  classical  analyses  have  been  largely  confined  to  networks  of  simple  shape  or 
of  geometry  that  lends  itself  to  analytical  or  semi-analytical  methods  of  solution. 

The  highly  complex  geometries  used  in  microwave  circuits  today  makes  it  necessary  to  use  numerical 
methods  for  analysis.  Multi-port  microwave  junctions  are  solved  numerically  in  the  literature  by  using 
one  of  two  approaches:  the  eigensolution  method3,3  and  the  deterministic  method7,3’3.  In  the  eigensolu- 
tion  method,  either  the  finite  element  method  or  the  finite  difference  method  is  used  to  compute  the 
eigenvalues  and  the  eigenvectors  of  the  normal  modes  of  the  junction,  and  then  circuit  theory  is  used  to 
determine  the  circuit  parameters10.  In  the  deterministic  method,  the  field  solution  is  computed  at  a  single 
specified  frequency  and  scattering  parameters  computed  at  that  frequency  only;  the  entire  process  must  be 
repeated  to  determine  the  solution  at  other  frequencies. 

While  the  eigensolution  method  is  mathematically  elegant,  it  has  the  disadvantage  of  requiring  the 
solution  of  large  matrix  eigenvalue  equations.  Since  the  solution  of  matrix  eigenvalue  problems  is  expen¬ 
sive,  recent  work  has  focused  on  the  deterministic  approach  that  requires  the  solution  of  deterministic 
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matrix  equations  only7’*.  Recently,  Webb*  has  used  finite  element  method  for  the  analysis  of  //-plane 
rectangular  waveguide  problems.  In  his  procedure,  a  set  of  boundary  value  problems  is  solved  in  order  to 
get  the  field  solution  at  a  single  frequency  point.  Two  similar  procedures  have  been  developed  by  Kon- 
shiba.  In  Reference  7,  the  boundary  element  is  combined  with  modal  analysis  to  solve  waveguide  discon¬ 
tinuity  problems  and  the  finite  element  method  is  combined  with  modal  analysis  to  solve  for  fields  in  an 
//-plane  waveguide  circulator  in  Reference  9.  Unfortunately,  however,  both  of  these  procedures  result  in 
non-symmetric  matrix  equations  that  are  expensive  to  solve. 

In  this  paper,  we  introduce  a  new,  highly  efficient  procedure  for  modeling  N- port  waveguide  junctions. 
The  basis  of  the  procedure  is  the  transfinite  element  method11,13  in  which  modal  basis  functions  are  com¬ 
bined  with  finite  element  basis  functions  to  provide  solutions  for  open  boundary  problems.  This  proce¬ 
dure  results  in  symmetric  sparse  matrix  equations  that  can  be  solved  very  efficiently  by  using  the  pre¬ 
conditioned  bi-conjugate  gradient  algorithm.  Further,  we  develop  a  spectral  response  estimation  proce¬ 
dure  by  which  solutions  at  a  few  adaptively  selected  frequencies  are  used  to  generate  the  full  solution  in 
the  frequency  range  of  interest. 

Numerical  results  are  given  for  a  T-junction,  a  screen  filter  containing  E-plane  metal  inserts  and  a 
dielectric  filter  to  show  the  validity  of  the  present  procedure.  A  comparison  of  the  computation  times 
required  by  the  adaptive  procedure  and  by  the  direct  deterministic  procedure  for  the  T-junction  problem 
is  also  presented. 

FORMULATION 

THE  TRANSFINITE  ELEMENT  METHOD 

The  structure  to  be  analysed  consists  of  a  cavity  coupled  with  N  rectangular  waveguides.  The  shape 
of  the  cavity  and  the  dimensions  of  the  waveguides  are  arbitrary,  but  the  overall  structure  must  be 
uniform  so  that  the  problem  can  be  approximated  by  two-dimensional  analysis.  To  simplify  the  formula¬ 
tion,  we  assume  that  the  junction  is  in  the  //-plane:  problems  involving  E-plane  junctions  can  be  treated 
in  much  the  same  way.  Microwave  planar  circuit  problems  differ  only  in  that  the  electric  field  is  taken  to 
be  a  constant  perpendicular  to  the  plane  of  the  circuit  and  that  regions  of  different  dielectric  constant  are 
seldom  used. 

Consider  exciting  port  1  by  the  dominant  TEJ0  mode.  The  field  over  the  eavity  region  f2.  and  the 
port  regions  extending  to  infinity  must  satisfy  the  Helmholts  equation: 
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V*E+  &trE  —  0 
w2  <0  ^0 

(1) 

where  u  it  the  angular  frequency  of  the  excitation,  (g  and  are  the  permittivity  and  the  permeability  of 
the  free  (pace,  respectively,  and  tf  it  the  relative  dielectric  constant  of  the  material. 

In  the  transfmite  element  method11’13,  the  problem  region  is  divided  into  two  parts.  An  interior 
region  ft.  of  finite  extent  and  an  exterior  region  that  is  homogeneous  and  unbounded.  Within  f?(.  finite 
element  basis  functions  are  used  to  approximate  the  field;  in  analytical  solutions  of  the  Helmholtz 
equation  provide  a  basis  set  for  the  field.  Employing  both  sets  of  basis  functions  in  a  variational  proce¬ 
dure  and  requiring  continuity  along  the  boundary  between  and  1?(  gives  a  symmetric  matrix  equation 
that  is  solved  for  the  field. 

The  application  of  the  transfmite  element  method  to  multi-port  microwave  circuits  is  presented  in 
Reference  12.  In  this  procedure,  Expiation  (1)  is  converted  into  the  matrix  equation 


<lS)-»*|Tf+l*l)  *  --T+**7 


(2) 

A  A  A 

where  [5]  and  [7]  are  complex  symmetric  matrices,  pi]  is  a  diagonal  matrix,  tP  is  the  solution  vector,  k  is 
the  wavenumber,  and  /  and  g  are  known  vectors. 


With  respect  to  spectral  modeling,  we  note  that  the  matrices  [5]  and  (T)  and  the  vectors  /  and  g  are 

frequency  independent  and  can  be  computed  mice  and  stored  for  any  problem  geometry.  Only  the 

wavenumber  k  and  the  matrix  (y]  depend  on  frequency.  However,  since  [y]  is  a  diagonal  matrix  with  only 

N  times  M  non- zero  entries,  where  N  is  the  number  of  ports  in  the  circuit  and  M  is  the  number  of  basis 

A 

functions  used  to  approximate  the  fields  in  each  port,  it  requires  very  little  work  to  compute  [y]. 


Figure  1  provides  the  transfmite  element  solution  of  a  microwave  T-junction  at  S  GHz. 
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SPECTRUM  MODELING 

Id  designing  microwave  components,  the  frequency  response  over  a  given  frequency  range  is  often  re¬ 
quired.  With  the  deterministic  approach  in  which  Equation  (1)  is  solved  at  a  given  frequency,  Equation 
(2)  must  be  solved  L  times  to  get  L  points  on  the  spectrum.  For  problems  where  the  response  changes 
▼ery  quickly,  the  number  of  points  L  used  to  generate  the  spectrum  must  be  large  in  order  to  get  satis¬ 
factory  answers.  The  adaptive  scheme  proposed  here  for  modeling  the  spectral  response  is  best  explained 
by  referring  to  Figure  2.  Instead  of  solving  Equation  (2)  L  times,  we  solve  the  matrix  equation  at  a  few 
adaptively  selected  optimal  frequencies  and  then  use  these  solutions  as  basis  functions  to  generate  the 
entire  spectral  response. 

A  preliminary  illustration  of  the  procedure  is  as  follows:  In  Figure  2(a),  the  three-port  junction  of 
Figure  1  is  solved  by  using  the  transfinite  element  method  at  the  two  limiting  frequencies  indicated  by  the 
dots  near  the  horisontal  axis.  These  two  solutions  are  then  used  as  basis  functions  to  generate  a  crude 
spectra]  response  curve  throughout  the  region  of  interest.  This  is  plotted  as  solid  lines  in  Figure  2(a). 
Next,  we  compute  the  error  in  the  solution  throughout  the  frequency  range  by  substituting  the  crude 
solution  values  into  the  governing  equations  for  the  system  and  by  evaluating  the  residual.  As  explained 
below,  this  may  be  done  very  efficiently.  We  then  solve  the  system  once  again  using  the  transfinite  ele¬ 
ment  algorithm  at  the  frequency  that  gave  the  maximum  residual  on  the  last  pass.  The  new  dot  in 
Figure  2(b)  shows  the  location  of  this  solution  as  well  as  the  new  spectral  response  curve  computed  by 
using  the  three  transfinite  element  solutions  as  basis  functions.  This  process  is  repeated  for  six  iterations 
in  Figure  2  until  the  error  in  the  entire  spectral  response  curve  is  within  acceptable  limits.  As  is  evident 
from  Figure  2(f),  the  procedure  converges  to  the  solution  given  in  Reference  4,  but  as  shown  below,  in 
much  less  computer  time. 

THE  ADAPTIVE  ALGORITHM 

To  develop  the  adaptive  spectral  response  modeling  algorithm,  first  approximate  the  solution  IP  of 
Equation  (2)  at  an  arbitrary  frequency  u  by  a  linear  combination  of  the  basis  functions  Ip . 


*  m  -  23  *  i 

1 


(3) 


where  n  is  the  number  of  basis  functions  and  a.  are  unknown  coefficients.  The  basis  functions  |P(.  are 
taken  to  be  the  solutions  of  (2)  at  n  specified  frequencies  Substituting  Equation  (3)  into  (2)  gives 
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4  A  a  ^ 

<[S|-*’l7]  +  M)*£ •/«)  -7  +*2 7 


»— 1 


(<) 

A  A  ***  ^ 

Since  [S],  [T],  / ,  and  g  are  frequency  independent,  we  can  rewrite  Equation  (4)  as 


£  «/w)(  +  (*]•  tfi  )  -  7  -**7 


(5) 


where 


7.  -  (S]*T 
7,.-  { 7j» 


(8) 

Now  applying  the  method  of  least-squares  to  Equation  (5)  leads  to  the  following  matrix  equation  for  the 
coefficients  a. 


where 


(7) 
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4.V*,+ 

£ 

4>7 

£%  7  +  £*•'? 

£•1 

Jt.-lTT]  /  -**  *,M  * 


Notice  that  in  (f),  only  the  matrix  [D]  and  the  rector  4  depend  on  the  frequency  <•/.  The  computation 

A 

time  for  evaluating  [D]  and  d  for  each  frequency  ia  negligible  because  [7]  is  a  diagonal  matrix.  Matrix 
Equation  (7)  is  therefore  trivial  to  set  up  and  it  is  inexpensive  to  solve  since  it  is  only  of  order  n. 

The  residua]  associated  with  the  frequency  u>  after  solving  for  the  coefficients  a^w)  is  evaluated  as 
follows 


A.  ~ 


7  -  7-*,7-r*<M<S'-*,#r+M  *,) 


The  norm  of  the  residual  is  given  by 
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K«)l  —  (  p**  7  ) i//2 

(10) 

The  entire  adaptive  solution  process  is  presented  in  the  flowchart  in  Figure  3.  First,  input  the  desired 
frequency  limits  uf,  and  L,  the  number  of  equal  divisions  in  the  frequency  range  to  be  used.  Next  solve 

(2)  at  the  two  limiting  frequencies  to  generate  the  basis  functions  and  the  auxiliary  functions  4>  , 

*v  ev 

+  2’  *v  & 2'  ^rom  functions  and  the  auxiliary  functions,  compute  the  components  of  matrices 

#>S* 

[A],  \B\,  \C\  and  of  vectors  F,  C,  H.  Then  generate  L  approximate  solutions  at  the  intermediate  fre- 
quencies  and  compute  the  residue  at  each  frequency  to  indicate  the  corresponding  error.  At  the  maximum 
residual,  solve  Equation  (2)  again  and  add  the  new  solution  to  the  basis  set  for  generating  the  spectrum. 
Repeat  the  process  until  the  maximum  residual  is  smaller  than  the  prescribed  error  tolerance  q. 

NUMERICAL  RESULTS 

Figures  4-5  provide  examples  of  spectrum  response  obtained  by  using  the  new  adaptive  spectrum 
modeling  procedure.  The  number  of  modes  used  in  each  port  to  generate  the  results  is  3  for  all  of  the 
problems  shown.  Solutions  of  Equation  (2)  are  obtained  in  the  computer  program  by  using  the  precon¬ 
ditioned  bi-conjugate  gradient  method  and  have  relative  residual  norms  smaller  than  I.0e-5.  The  error 
tolerance  q  employed  is  1.0e-4. 

Figure  2  shows  the  convergence  of  the  procedure  for  modeling  the  T-junction  in  Figure  1  throughout 
the  frequency  range  of  dominant  mode  propagation.  The  number  of  frequency  points  L  used  was  50,  and 
the  frequencies  that  are  solved  for  and  used  as  the  basis  functions  in  the  procedure  are  indicated  by  dots 
on  the  axis.  As  shown  in  Figure  2,  the  procedure  converges  when  n  6.  Comparing  the  final  adaptivity 
produced  spectral  response  with  the  boundary  element  deterministic  solution  of  Reference  7  shows  excel¬ 
lent  agreement. 

The  analysis  presented  here  is  not  limited  to  dominant  mode  propagation.  Figure#  shows  the  spectral 
response  computed  by  the  adaptive  spectrum  modeling  procedure  for  the  £-plane  metal  insert  filter  shown 
in  Figure  4(f)  from  30  GHz  to  60  GHz.  The  number  of  frequency  points  used  was  100  for  the  plots  in 
Figure  4(a)  -  4(e),  and  the  procedure  terminates  in  six  iterations.  A  comparison  of  the  final  spectrum 
response  with  that  obtained  by  the  field  expansion  calculation  in  Reference  13  is  given  in  Figure  4(e). 


177 


Figure  5(a)  presents  the  adaptive  solution  of  a  waveguide  dielectric  filter  modeled  by  Koshiba  and 
Sutuki7.  Here  we  used  L  ■«  100  spectral  points;  the  total  number  of  solved  frequency  points  was  only  5. 

Figure  6  shows  a  comparison  of  the  computer  time  required  by  the  new  procedure  and  that  of  the 
deterministic  approach.  The  time  for  the  deterministic  approach  is  based  on  the  total  time  needed  to 
solve  Equation  (2).  The  times  reported  are  for  the  ^junction  problem  using  a  DEC  VAX  11/780  com¬ 
puter;  the  matrix  site  was  500  by  500. 

CONCLUSIONS 

A  very  efficient  procedure  for  determining  the  spectral  response  of  microwave  circuits  has  been 
developed.  The  procedure  may  be  applied  to  waveguide  junctions  involving  either  E-plane  or  //-plane 
discontinuities  or  to  microwave  planar  circuits.  The  spectral  response  evaluation  procedure  employs  the 
transfinite  element  method  to  solve  for  the  field  at  a  few  adaptively  selected  frequencies  and  then  con¬ 
structs  the  solution  at  any  frequency  by  using  the  computed  solutions  as  basis  functions.  In  typical 

blems,  only  five  or  six  transfinite  element  solutions  are  required  to  converge  to  the  full  spectral 
response  evaluated  at  100  points  throughout  the  frequency  range  of  interest. 

In  the  past,  there  were  two  basic  alternatives  to  computing  the  spectral  response  of  microwave  circuits: 
One  could  employ  the  eigenaolution  approach  that  required  expensive  eigenvalue  problems  to  be  solved 
but  gave  solutions  at  in-between  frequencies  very  economically,  or  one  could  employ  the  existing  deter¬ 
ministic  approach  that  provided  solutions  at  a  specified  frequency  relatively  efficiently  but  had  to  be 
reapplied  at  every  frequency  of  interest.  The  new  spectral  response  modeling  procedure  combines  the 
advantages  of  both  approaches  and  gives  the  full  spectral  response  in  orders  of  magnitude  less  computing 
time  than  the  alternatives. 
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Figure  3: 


Power  Reflection  and  Trans  minion  Coefficients  of  tbe  T-  Junction  in  Figure  1. 
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Figur*  S: 


Flow  Chart  of  the  Adaptive  Spectrum  Modeling  Procedure. 


Figure  4: 


Power  Trusmianon  Coefficient  of  ea  S-Plase  Iatefrated  Screes  Filter. 
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